Zdrojový kód k článku:

• GitHub (zdrojový kód): https://github.com/Semtexcz/python-modeling-reality/blob/main/notebooks/vector_fields.ipynb
• Google Colab (spustit online): https://colab.research.google.com/drive/1r8G2HgYz01bzO_HKE1h8rg7QfiiMa4Wl?usp=sharing

Na vysoké škole jsem se s matematikou docela trápil. Přišel jsem z chemické průmyslovky, kde se na matematiku nekladl takový důraz jako na gymnáziích. První rok studia jsem proto ve velké míře doháněl matematické základy. Počítal jsem velké množství příkladů a často jsem se matematiku učil spíše mechanicky než skutečným porozuměním.

Hlavní problém byl v tom, že mi tehdy mnoho věcí nedávalo smysl. Často jsem si kladl otázku, kterou si klade spousta studentů:
„K čemu mi to vlastně bude?“

Musím říct, že způsob výuky tomu někdy příliš nepomáhal. Výklad často probíhal velmi úsporně a formálně: definice, věta, důkaz, jeden ukázkový příklad a pokračujeme dál, protože sylabus je dlouhý a času málo. Tento styl výuky je sice efektivní z hlediska předání formální struktury matematiky, ale neodpovídá tomu, jak matematika ve skutečnosti historicky vznikala.

Historický proces objevování matematiky je totiž téměř opačný.

Většinou začíná nějakým problémem. Člověk, skupina lidí nebo celé lidstvo narazí na konkrétní otázku. Často jde o velmi specifický problém – matematik by řekl speciální případ. Lidé začnou problém zkoumat a postupně si všímají určitých pravidelností nebo vzorů.

Například si někdo může všimnout, že obvod kruhu roste úměrně s jeho průměrem. Pokud změří několik kruhů, zjistí, že poměr obvodu a průměru je vždy přibližně stejný. Nejprve vyjde číslo okolo 3. Později přesnější měření dávají hodnoty 3,1; 3,14; 3,141; 3,1415 a tak dále.

Teprve ve chvíli, kdy je tento vzor dostatečně zřejmý, vznikne matematické zobecnění:

$O = \pi d$

nebo

$O = 2\pi r$

Matematický vzorec tak vlastně vzniká ze vzoru, který jsme v realitě pozorovali. Slovo „vzorec“ tuto souvislost krásně připomíná.

Matematika tedy historicky nevznikala tak, že by někdo nejprve napsal definici, poté větu a důkaz. Většinou to bylo naopak: nejdřív přišlo pozorování a hledání pravidelností a až následně jejich přesné matematické vyjádření.

A právě tento způsob uvažování mi začal dávat smysl až mnohem později, když jsem matematiku začal vnímat jako jazyk pro popis vzorů ve světě kolem nás. V tomto článku bych chtěl jeden takový jazyk představit: matematiku vektorových polí, která stojí za popisem proudění tekutin, elektromagnetismu, difuze nebo přenosu tepla.

Ironii osudu je, že jsem na tomto tématu vyletěl poprvé a naposledy z písemné zkoušky. Tehdy jsem ještě nechápal, k čemu se to používá, aplikace přišla později.

Co je pole?

Slovo pole má v češtině několik významů. Většina lidí si pravděpodobně představí pole, na kterém se pěstují plodiny. Pokud programujete, možná vás napadne také pole jako datová struktura, která slouží k ukládání sekvence hodnot v paměti počítače.

V tomto článku ale nemám na mysli ani jedno z těchto polí. Budeme se bavit o matematickém poli. Než se k němu dostaneme formálněji, podívejme se na obrázek, který jeho pochopení výrazně usnadní.

Na obrázku vidíme velkou část mapy světa. Zobrazuje dvě informace, které jsou pro nás důležité.

První informací je barva jednotlivých oblastí. Barva zde reprezentuje teplotu vzduchu v daném místě. Každému bodu na mapě je tedy přiřazeno jedno číslo – teplota. Takový typ matematického objektu nazýváme skalární pole.

Formálně řečeno je skalární pole funkce, která každému bodu v prostoru přiřazuje jednu hodnotu, tzv. skalár. Typickými příklady skalárních polí jsou například:

• teplota v atmosféře,
• tlak v kapalině,
• koncentrace látky v chemickém reaktoru,
• nebo nadmořská výška na mapě terénu.

Druhou informaci na obrázku představují bílé čáry, které vypadají jako jemné šmouhy. Ty znázorňují směr proudění větru. V každém bodě prostoru zde tedy není jen jedna hodnota, ale informace o směru a velikosti pohybu vzduchu.

Takový typ pole nazýváme vektorové pole.

Formálně je vektorové pole funkce, která každému bodu v prostoru přiřazuje vektor. Vektor je veličina, která má nejen velikost, ale také směr. Proto se hodí k popisu jevů, kde záleží na orientaci v prostoru.

Typickými příklady vektorových polí jsou například:

• rychlost proudění vzduchu (vítr),
• rychlost proudění vody v řece,
• elektrické pole kolem nabité částice,
• magnetické pole kolem vodiče s proudem.

Skalární pole nám tedy říká kolik, zatímco vektorové pole nám říká kolik a kam. Právě tato jednoduchá myšlenka stojí za velkou částí matematiky, která popisuje proudění tekutin, elektromagnetismus nebo difuzi látek.

Pokud chcete vidět tu mapu interaktivně, koukněte sem: https://www.ventusky.com/

Pole v Pythonu

Můj blog se snažím psát tak, aby byl co nejvíce lidský, konkrétní a praktický. Nechci jen vysvětlovat matematické pojmy abstraktně, ale ukázat také, jak s nimi můžeme pracovat v praxi. Proto v článku najdete i jednoduché ukázky kódu v Pythonu, které si můžete sami vyzkoušet.

Začněme tím nejjednodušším typem pole, se kterým se v matematice i fyzice setkáme – skalárním polem. Na krátkém příkladu si ukážeme, jak takové pole vytvořit a vizualizovat.

Skalární pole

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-2,2,100)
y = np.linspace(-2,2,100)

X, Y = np.meshgrid(x,y)

Z = X**2 + Y**2

plt.contourf(X,Y,Z,20)
plt.colorbar()
plt.title("Scalar field")
plt.show()

Pojďme si postupně vysvětlit, co jednotlivé části kódu dělají.

Nejprve importujeme dvě knihovny:

import numpy as np
import matplotlib.pyplot as plt

Knihovna NumPy slouží k numerickým výpočtům a práci s poli čísel.
Matplotlib používáme pro vizualizaci dat.

Vytvoření souřadnic

x = np.linspace(-2,2,100)
y = np.linspace(-2,2,100)

Funkce linspace vytvoří 100 rovnoměrně rozložených bodů mezi hodnotami -2 a 2.
Tyto body budou představovat souřadnice v rovině.

Můžeme si je představit jako jemnou síť bodů na mapě.

Vytvoření mřížky bodů

X, Y = np.meshgrid(x,y)

Funkce meshgrid vytvoří dvourozměrnou mřížku bodů ze souřadnic x a y.

Každý bod mřížky má tedy souřadnice:

$(x,y)$

Tato mřížka reprezentuje prostor, ve kterém budeme skalární pole počítat.

Definice skalárního pole

Z = X**2 + Y**2

Zde definujeme funkci

$ f(x,y) = x^2 + y^2 $

Tato funkce každému bodu roviny přiřazuje jedno číslo.

Například:

To je přesně definice skalárního pole.

Vykreslení pole

plt.contourf(X,Y,Z,20)

Funkce contourf vykreslí vyplněné vrstevnice (contour plot).
Číslo 20 určuje počet barevných úrovní.

Výsledek připomíná topografickou mapu, kde různé barvy reprezentují různé hodnoty funkce.

Barevná legenda

plt.colorbar()

Přidá barevnou škálu, která ukazuje, jaké hodnoty odpovídají jednotlivým barvám.

Titulek grafu

plt.title("Scalar field")

Přidá název grafu.

Zobrazení grafu

plt.show()

Nakonec graf vykreslí na obrazovku.

Co vlastně vidíme?

Výsledný obrázek ukazuje skalární pole funkce

$ f(x,y) = x^2 + y^2 $

Nejmenší hodnota je uprostřed (v bodě (0,0)).
Čím dál od středu se posouváme, tím hodnota funkce roste.

Geometricky jde o paraboloid – tvar připomínající misku.

Teď se posuňme o krok dál a podívejme se na vektorová pole. Na rozdíl od skalárního pole, které každému bodu přiřazuje jednu hodnotu, vektorové pole přiřazuje vektor – tedy veličinu, která má velikost i směr.

Abychom si tuto myšlenku lépe představili, ukážeme si jednoduchý příklad v Pythonu. Pomocí několika řádků kódu můžeme vektorové pole nejen spočítat, ale také vizualizovat jako síť šipek, které ukazují směr a velikost vektorů v jednotlivých bodech prostoru.

Vektorové pole

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-2,2,20)
y = np.linspace(-2,2,20)

X, Y = np.meshgrid(x,y)

U = -Y
V = X

plt.quiver(X,Y,U,V)
plt.title("Vector field")
plt.show()

Pojďme si vysvětlit jednotlivé části.

Import knihoven

Nejprve importujeme dvě knihovny, stejné jako u skalárního pole.

import numpy as np
import matplotlib.pyplot as plt

Vytvoření souřadnic

x = np.linspace(-2,2,20)
y = np.linspace(-2,2,20)

Funkce linspace vytvoří 20 rovnoměrně rozložených bodů mezi hodnotami -2 a 2.

Tyto body představují souřadnice v rovině. V podstatě si připravujeme mřížku bodů, ve kterých budeme vektorové pole počítat.

Vytvoření mřížky

X, Y = np.meshgrid(x,y)

Funkce meshgrid vytvoří dvourozměrnou síť bodů.

Každý bod má souřadnice

$(x,y)$

Tato mřížka reprezentuje prostor, ve kterém budeme vektorové pole definovat.

Definice vektorového pole

Teď definujeme samotné pole.

U = -Y
V = X

Každý vektor v poli má dvě složky:

$F(x,y) = (U,V)$

V našem případě tedy

$F(x,y) = (-y, x)$

To znamená:

• vodorovná složka vektoru je (-y)
• svislá složka je (x)

Jak pole vypadá

Toto pole má zajímavou vlastnost: vektory rotují kolem středu.

Například:

Vektory tedy tvoří rotační pole – podobné například víru v kapalině.

Vykreslení pole

Pole vykreslíme pomocí funkce quiver.

plt.quiver(X,Y,U,V)

Funkce quiver zobrazí vektorové pole jako síť šipek.

Každá šipka:

• začíná v bodě ((x,y))
• směřuje podle vektoru ((U,V))

Titulek a zobrazení

Nakonec přidáme název grafu a vykreslíme ho.

plt.title("Vector field")
plt.show()

V předchozím příkladu jsme použili funkci quiver, která zobrazuje vektorové pole jako mřížku šipek. Každá šipka ukazuje směr a velikost vektoru v daném bodě prostoru. Pro pochopení struktury pole je to velmi užitečné, protože vidíme jednotlivé vektory přímo.

Existuje ale ještě jiný způsob vizualizace vektorového pole, který je často ještě intuitivnější. Místo jednotlivých šipek můžeme zobrazit tzv. proudnice (anglicky streamlines). Ty ukazují, jak by se v daném poli pohybovala částice, kdyby se v něm nechala unášet. V Pythonu k tomu slouží funkce streamplot z knihovny Matplotlib. Výsledný obrázek pak připomíná například proudění kapaliny, větru nebo magnetické siločáry. Níže je ukázka kódu a výsledného obrázku.

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-3,3,100)
y = np.linspace(-3,3,100)

X, Y = np.meshgrid(x,y)

U = -Y
V = X

speed = np.sqrt(U**2 + V**2)

plt.streamplot(X, Y, U, V, color=speed, cmap="viridis")
plt.colorbar(label="Speed")
plt.title("Vector field")
plt.show()

Jak můžeme pole popisovat?

Zavedli jsme si pojem pole a ukázali jsme si dva základní typy: skalární pole a vektorové pole. Viděli jsme také, jak taková pole můžeme jednoduše vizualizovat pomocí Pythonu.

A co dál?

Pokud chceme pomocí matematiky modelovat reálné jevy, nestačí pole jen definovat. Musíme také umět popsat, jak se pole v prostoru mění, kde vznikají jejich maxima a minima, kde se koncentrují zdroje nebo kde dochází k proudění.

Abychom si tyto pojmy představili intuitivně, vraťme se na chvíli k příkladu počasí. Z předpovědi počasí pravděpodobně znáte pojmy jako tlaková výše, tlaková níže nebo teplá a studená fronta. Meteorologové jimi popisují právě strukturu různých polí v atmosféře – například pole tlaku, teploty nebo rychlosti větru.

Tlaková výše a níže například odpovídají místům, kde má pole tlaku lokální maxima nebo minima, zatímco fronty označují oblasti, kde se hodnoty polí rychle mění v prostoru.

V matematice existují nástroje, které nám umožňují tyto změny přesně popsat. Nejčastěji se používají tři základní operace nad poli: gradient, divergence a rotace. V následujících částech si ukážeme, co tyto pojmy znamenají a jak souvisejí s fyzikálními jevy, které známe z každodenního života.

Gradient

Představte si následující situaci. Jedete na kole. V prvním případě jedete po cyklostezce v Nizozemsku podél nějakého kanálu. Krajina je téměř dokonale rovná a změny nadmořské výšky jsou minimální. Ve druhém případě jedete v Italských Alpách a snažíte se vyšlapat legendární horský průsmyk Passo Stelvio.

Na kterém z těchto míst bude nadmořská výška růst rychleji?

Otázka je samozřejmě řečnická. V Alpách budete stoupat mnohem prudčeji než na rovinaté nizozemské cyklostezce.

V matematice a fyzice bychom tuto situaci popsali pomocí gradientu. Gradient nám říká, jak rychle se nějaká veličina mění v prostoru a kterým směrem tato změna probíhá nejrychleji.

Pokud bychom nadmořskou výšku popsali jako skalární pole $h(x,y)$, pak by gradient tohoto pole ukazoval směr, ve kterém výška roste nejrychleji. Velikost gradientu pak říká, jak prudké toto stoupání je.

V našem příkladu tedy můžeme říct, že gradient nadmořské výšky ve směru jízdy je v Italských Alpách mnohem větší než na rovinaté cyklostezce v Nizozemsku. Právě proto se v Alpách nadřete mnohem více než na rovině.

Příklad s kolem nám dává dobrou intuici. Vidíme, že nadmořská výška se může v různých místech měnit různě rychle. Na rovině téměř vůbec, zatímco v horách může stoupat velmi prudce.

Teď si představte jinou situaci. Stojíte na svahu na lyžích a vaším cílem je dosáhnout co největší rychlosti. Musíte tedy zvolit takový směr, ve kterém je kopec nejprudší. Zkušení lyžaři by řekli, že se máte vydat po spádnici.

A co je vlastně spádnice?

Spádnice je směr, ve kterém nadmořská výška klesá (nebo roste, pokud jdeme opačným směrem) nejrychleji. Jinými slovy, je to směr největší změny nadmořské výšky v daném bodě svahu.

Právě tento směr popisuje v matematice gradient.

Pokud nadmořskou výšku zapíšeme jako funkci dvou proměnných

$$h(x,y)$$

pak gradient této funkce definujeme jako

$$\nabla h = \left(\frac{\partial h}{\partial x},\frac{\partial h}{\partial y} \right)$$

Tento zápis říká dvě věci zároveň. První složka popisuje, jak rychle se nadmořská výška mění ve směru osy $x$, druhá složka říká, jak rychle se mění ve směru osy $y$. Když tyto dvě informace spojíme dohromady, dostaneme vektor, který ukazuje směr nejrychlejšího růstu výšky – tedy přesně to, co jsme na svahu intuitivně nazvali spádnicí.

Gradient tedy převádí skalární pole (například pole nadmořské výšky) na vektorové pole, které v každém bodě ukazuje směr největší změny. Právě proto se gradient objevuje v mnoha fyzikálních zákonech – například při popisu proudění tepla, difuze látek nebo elektrického pole.

Python ukázka gradientu

Vrátíme se ke skalárnímu poli, které jsem již v článku použili $f(x, y) = x^2 + y^2$ a ukážeme si na něm jak spočítat gradient a zobrazit ho jako vektorové pole.

import numpy as np
import matplotlib.pyplot as plt

# vytvoření souřadnic
x = np.linspace(-2, 2, 40)
y = np.linspace(-2, 2, 40)

X, Y = np.meshgrid(x, y)

# definice skalárního pole
Z = X**2 + Y**2

# výpočet gradientu
dZdy, dZdx = np.gradient(Z)

# vizualizace
plt.contourf(X, Y, Z, 20, cmap="viridis")
plt.colorbar(label="Value")

plt.quiver(X, Y, dZdx, dZdy, color="white")

plt.title("Gradient of scalar field")
plt.show()

Co tento kód dělá

Nejprve vytvoříme mřížku bodů v rovině pomocí funkcí linspace a meshgrid.

X, Y = np.meshgrid(x, y)

Každý bod této mřížky má souřadnice ((x,y)).

Definice skalárního pole

Pak definujeme funkci

$$f(x,y) = x^2 + y^2$$

v Pythonu:

Z = X**2 + Y**2

Toto pole můžeme interpretovat například jako nadmořskou výšku kopce.

Výpočet gradientu

Gradient spočítáme pomocí funkce np.gradient.

dZdy, dZdx = np.gradient(Z)

Tato funkce vrací dvě pole:

• změnu ve směru (x)
• změnu ve směru (y)

tedy právě složky gradientu.

Vizualizace

Nejprve vykreslíme skalární pole jako barevnou mapu.

plt.contourf(X, Y, Z, 20)

Potom přes něj vykreslíme gradient pomocí šipek.

plt.quiver(X, Y, dZdx, dZdy)

Každá šipka ukazuje směr nejrychlejšího růstu funkce.

Co vidíme na obrázku

Výsledný obrázek ukazuje zajímavou vlastnost gradientu:

• všechny šipky míří od středu ven
• uprostřed je minimum funkce

To odpovídá intuici – pokud stojíme v nejnižším bodě „údolí“, nejrychleji začneme stoupat směrem ven od středu.

Výsledný obrázek není vizuálně úplně působivý, proto ho můžeme drobně upravit.

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)

X, Y = np.meshgrid(x, y)

Z = X**2 + Y**2

dZdy, dZdx = np.gradient(Z)

plt.contourf(X, Y, Z, 30, cmap="terrain")
plt.colorbar(label="Elevation")

step = 6

plt.quiver(
    X[::step, ::step],
    Y[::step, ::step],
    dZdx[::step, ::step],
    dZdy[::step, ::step],
    color="black",
    scale=None,          # vypne automatické škálování
    width=0.004
)

plt.title("Gradient of scalar field (elevation)")
plt.xlabel("x")
plt.ylabel("y")

plt.show()

Divergence

Než budeme pokračovat podívejte se prosím na toto video.

Ve videu vidíme zapálenou dýmovnici a vedle ní čističku vzduchu, která nasává okolní vzduch a filtruje ho. Z hlediska analýzy polí je to velmi pěkný a názorný experiment.

Dýmovnice uvolňuje do prostoru velké množství prachových částic. Ty se začnou pohybovat vzduchem a v každém bodě prostoru mají určitou rychlost a směr pohybu. Jinými slovy, můžeme si představit, že v každém bodě existuje vektor rychlosti, který popisuje pohyb částic v daném místě. Vzniká tak vektorové pole proudění částic.

Jakmile zapneme čističku vzduchu, situace se změní. Částice jsou nasávány směrem k čističce, kde jsou zachyceny filtrem. V okolí dýmovnice se tedy částice rodí (jsou do prostoru uvolňovány), zatímco u čističky zanikají (jsou pohlceny filtrem).

Pokud to chceme popsat přesněji jazykem vektorové analýzy, pak nás nezajímá přímo vznik nebo zánik jednotlivých částic, ale to, jak se chová tok vektorů v prostoru. V některých místech se vektory rozbíhají do okolí, jinde se naopak sbíhají do jednoho bodu.

Právě tuto vlastnost pole popisuje divergence. Divergence v každém bodě říká, zda se tok vektorů v daném místě rozbíhá (zdroj), sbíhá (odtok, propad), nebo zda jím pouze protéká bez změny.

Slovo divergence pochází z latinského divergere. Je složené ze dvou částí: předpony dis- / di-, která znamená „od sebe“, a slovesa vergere, tedy „směřovat“ nebo „naklánět se“. Doslovný význam slova divergence je tedy rozbíhání se od sebe. To velmi dobře odpovídá matematickému významu: vektorové pole se v daném místě rozbíhá do okolí.

Opakem divergence je konvergence. Toto slovo pochází z latinského convergere, kde předpona con- znamená „spolu“ nebo „dohromady“. Konvergence tedy doslova znamená sbíhání se k sobě. V kontextu vektorových polí to znamená, že tok směřuje do jednoho bodu, kde se „shromažďuje“.

V matematice tyto dva pojmy souvisejí také se znaménkem divergence. Pokud je divergence kladná, vektory se z daného místa rozbíhají a můžeme mluvit o zdroji pole. Pokud je divergence záporná, vektory se naopak sbíhají do jednoho bodu – pole má v tomto místě pohlcovač (někdy se používá také termín propad). Pokud je divergence nulová, tok vektoru daným místem pouze protéká, aniž by zde vznikal nebo zanikal.

Navážu tak, aby přechod z intuice (dýmovnice, pohlcovač, rozbíhání toků) k definici byl přirozený a ne formální.

Intuitivní představu tedy máme: divergence nám říká, zda se tok vektorů v daném místě rozbíhá, sbíhá, nebo zda jím pouze protéká. Jak ale tuto myšlenku vyjádřit matematicky?

Představme si malý uzavřený objem prostoru, například malou krychli. Sledujeme, kolik toku do této krychle vstupuje a kolik z ní vystupuje. Pokud z krychle vystupuje více toku, než do ní vstupuje, znamená to, že se tok v tomto místě rozbíhá – máme kladnou divergenci. Pokud naopak více toku do krychle vstupuje, než z ní vystupuje, tok se zde sbíhá – divergence je záporná.

Divergence tedy v podstatě měří rozdíl mezi tokem, který z malého objemu vystupuje, a tokem, který do něj vstupuje.

Matematicky tuto myšlenku zapisujeme pomocí parciálních derivací. Pokud máme vektorové pole

$$ \mathbf{F}(x,y,z) = (F_x(x,y,z), F_y(x,y,z), F_z(x,y,z))$$

pak jeho divergence je definována jako

$$ \nabla \cdot \mathbf{F} = \frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial F_z}{\partial z}.$$

Každý člen tohoto výrazu říká, jak se mění příslušná složka vektorového pole ve svém směru. Když tyto změny sečteme, dostaneme informaci o tom, zda se pole v daném místě lokálně rozbíhá, nebo sbíhá.

Tato definice má jednu velmi pěknou interpretaci: divergence říká, kolik toku z malého objemu prostoru „odtéká“ ven. Proto se s ní setkáváme v mnoha fyzikálních zákonech – například při popisu proudění tekutin, difuze látek nebo elektrického pole.

Python ukázka divergence

V hodinách fyziky na základní nebo střední škole jste se určitě setkali s elektrickým polem. Jedno z nejjednodušších polí je pole, kde jsou dva náboje. Jeden kladný a druhý záporný. Mezi těmito náboji vzniká elektrické pole, což je velmi hezký příklad pro vizualizaci v Pythonu.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import TwoSlopeNorm

# mřížka prostoru
x = np.linspace(-2, 2, 200)
y = np.linspace(-2, 2, 200)
X, Y = np.meshgrid(x, y)

# pozice nábojů
q1 = 1
q2 = -1
x1, y1 = -1, 0
x2, y2 = 1, 0

# vzdálenosti od nábojů
eps = 0.05
r1 = np.sqrt((X - x1)**2 + (Y - y1)**2 + eps**2)
r2 = np.sqrt((X - x2)**2 + (Y - y2)**2 + eps**2)

# elektrické pole
Ex = q1 * (X - x1) / r1**3 + q2 * (X - x2) / r2**3
Ey = q1 * (Y - y1) / r1**3 + q2 * (Y - y2) / r2**3

# divergence
dx = x[1] - x[0]
dy = y[1] - y[0]
dExdx = np.gradient(Ex, dx, axis=1)
dEydy = np.gradient(Ey, dy, axis=0)
div = dExdx + dEydy

# ořez extrémních hodnot pro lepší vizualizaci
limit = np.percentile(np.abs(div), 98)
div_clipped = np.clip(div, -limit, limit)

# barevná normalizace se středem v nule
norm = TwoSlopeNorm(vmin=-limit, vcenter=0, vmax=limit)

# vykreslení
plt.figure(figsize=(9, 7))

plt.contourf(
    X, Y, div_clipped,
    levels=40,
    cmap="RdBu_r",
    norm=norm
)
plt.colorbar(label="Divergence")

plt.streamplot(X, Y, Ex, Ey, color="black", density=1.2, linewidth=1)

plt.scatter(x1, y1, color="red", s=120, label="Positive charge")
plt.scatter(x2, y2, color="blue", s=120, label="Negative charge")

plt.title("Electric field and divergence")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

Vytvoření prostoru

Nejprve si vytvoříme dvourozměrnou mřížku bodů, ve které budeme elektrické pole počítat.

x = np.linspace(-2, 2, 200)
y = np.linspace(-2, 2, 200)
X, Y = np.meshgrid(x, y)

Funkce linspace vytvoří rovnoměrně rozložené body mezi -2 a 2. Funkce meshgrid z nich sestaví síť bodů v rovině. Každý bod této sítě odpovídá jednomu místu v prostoru, kde budeme pole vyhodnocovat.

Definice nábojů

Potom definujeme dva bodové náboje:

q1 = 1
q2 = -1
x1, y1 = -1, 0
x2, y2 = 1, 0

První náboj je kladný a leží vlevo, druhý je záporný a leží vpravo. Tím vytvoříme klasickou dvojici, u které elektrické pole vychází z kladného náboje a vstupuje do záporného.

Vzdálenosti od nábojů

Abychom mohli spočítat pole v každém bodě, potřebujeme znát vzdálenost daného bodu od každého náboje.

eps = 0.05
r1 = np.sqrt((X - x1)**2 + (Y - y1)**2 + eps**2)
r2 = np.sqrt((X - x2)**2 + (Y - y2)**2 + eps**2)

Zde používáme běžný vzorec pro eukleidovskou vzdálenost. Malá konstanta eps slouží k tomu, aby ve výpočtu nevznikla singularita přímo v místě náboje. V ideální fyzice by bylo pole bodového náboje v jeho středu nekonečné, ale při numerickém výpočtu potřebujeme tuto singularitu mírně zjemnit.

Výpočet elektrického pole

Nyní spočítáme složky elektrického pole v každém bodě.

Ex = q1 * (X - x1) / r1**3 + q2 * (X - x2) / r2**3
Ey = q1 * (Y - y1) / r1**3 + q2 * (Y - y2) / r2**3

Tento zápis odpovídá tomu, že elektrické pole bodového náboje míří radiálně od náboje a jeho velikost klesá s rostoucí vzdáleností. Výsledné pole dostaneme jako součet polí obou nábojů.

Pole tedy v každém bodě reprezentujeme dvojicí hodnot:

$$\mathbf{E}(x,y) = (E_x(x,y), E_y(x,y))$$

Výpočet divergence

Divergence vektorového pole ve dvou rozměrech je definována jako

$$\nabla \cdot \mathbf{E} = \frac{\partial E_x}{\partial x} + \frac{\partial E_y}{\partial y}.$$

V Pythonu ji spočítáme numericky pomocí funkce np.gradient.

dx = x[1] - x[0]
dy = y[1] - y[0]
dExdx = np.gradient(Ex, dx, axis=1)
dEydy = np.gradient(Ey, dy, axis=0)
div = dExdx + dEydy

Nejprve si spočítáme krok mřížky ve směru os $x$ a $y$. Potom určujeme parciální derivace obou složek elektrického pole. Jejich součtem získáme divergenci.

Proč je potřeba oříznout extrémní hodnoty

V okolí bodových nábojů jsou hodnoty divergence numericky velmi vysoké. Kdybychom je zobrazili bez úprav, barevná mapa by byla téměř nečitelná a většina zajímavých struktur by zanikla.

Proto extrémní hodnoty ořízneme:

limit = np.percentile(np.abs(div), 98)
div_clipped = np.clip(div, -limit, limit)

Tím si ponecháme většinu informace, ale zabráníme tomu, aby několik extrémních bodů zničilo celou škálu.

Barevná škála se středem v nule

Protože divergence může být kladná i záporná, chceme, aby nulová hodnota ležela přesně uprostřed barevné škály.

norm = TwoSlopeNorm(vmin=-limit, vcenter=0, vmax=limit)

To je důležité, protože:

• kladná divergence odpovídá zdroji pole,
• záporná divergence odpovídá pohlcovači,
• nula znamená, že pole místem jen protéká.

Vykreslení výsledku

Nejprve zobrazíme divergenci jako barevné pole.

plt.contourf(
    X, Y, div_clipped,
    levels=40,
    cmap="RdBu_r",
    norm=norm
)
plt.colorbar(label="Divergence")

Použili jsme barevnou mapu RdBu_r, protože je vhodná pro veličiny se zápornými i kladnými hodnotami. Jedna barva reprezentuje kladné hodnoty, druhá záporné.

Potom přes tepelnou mapu vykreslíme proudnice elektrického pole:

plt.streamplot(X, Y, Ex, Ey, color="black", density=1.2, linewidth=1)

Díky tomu na jednom obrázku vidíme zároveň:

• vektorové pole elektrického pole,
• skalární pole jeho divergence.

Nakonec vykreslíme samotné náboje:

plt.scatter(x1, y1, color="red", s=120, label="Positive charge")
plt.scatter(x2, y2, color="blue", s=120, label="Negative charge")

Jak výsledek interpretovat

Na výsledném obrázku uvidíme, že:

• v okolí kladného náboje je divergence kladná, protože odtud pole vychází,
• v okolí záporného náboje je divergence záporná, protože tam pole vstupuje,
• ve většině ostatního prostoru je divergence blízká nule.

To velmi dobře odpovídá fyzikální intuici i Gaussovu zákonu:

$$ \nabla \cdot \mathbf{E} = \frac{\rho}{\varepsilon_0} $$

Divergence elektrického pole je tedy nenulová právě tam, kde se nachází elektrický náboj.

Rotace

Na začátku jsme pracovali s mapou počasí. Když si zobrazíme skalární pole atmosférického tlaku na Zemi a zároveň vektorové pole větru, můžeme si všimnout zajímavého jevu, který zatím nemáme jak popsat.

Jinými slovy, vzduch se zde nepohybuje jen z jednoho místa do druhého. Celé proudění má rotační charakter – částice vzduchu opisují zakřivené dráhy a proudnice se kolem center tlaku stáčejí.

Právě tento jev — tedy lokální otáčení vektorového pole — popisuje další operátor vektorové analýzy, který se nazývá rotace (anglicky curl).

Rotace nám říká, zda má pole v daném místě tendenci vytvářet vír. Pokud bychom si představili malou lopatku nebo vrtulku umístěnou do proudícího pole, rotace nám říká, zda by se tato lopatka začala otáčet, a pokud ano, kolem jaké osy.

V meteorologii se rotace používá například při popisu cyklon a anticyklon, tedy velkých atmosférických vírů spojených s tlakovými nížemi a výšemi.

Matematika tak dokáže velmi přesně popsat jev, který známe z každodenní zkušenosti — vítr, který se stáčí do spirály kolem tlakových útvarů v atmosféře.

Divergence nám říká, zda se pole v daném místě rozbíhá nebo sbíhá. Viděli jsme to například u elektrického pole nebo proudění částic.

U větru na mapě počasí jsme ale pozorovali ještě jiný jev. Proudnice se kolem tlakových útvarů stáčejí do spirály. Jinými slovy, pole má vířivý charakter.

Jak ale takový vír matematicky změřit?

Představme si jednoduchý experiment. Do proudění položíme malé kolečko nebo lopatku, podobně jako malý vodní mlýnek. Pokud proudění působí na kolečko tak, že ho začne otáčet, znamená to, že pole má v daném místě rotaci.

Teď si představme velmi malý čtverec v prostoru. Sledujme, jak působí vektorové pole na jednotlivé jeho strany.

Pokud by rychlost pole byla na všech stranách stejná, čtverec by se jen posouval — ale neotáčel by se. Aby vznikla rotace, musí být situace trochu jiná.

Například:

• na horní straně čtverce proudění táhne doprava silněji,
• na dolní straně táhne doprava slaběji.

V takové situaci vznikne otáčivý moment a čtverec se začne otáčet.

Podobně může vzniknout rotace i v druhém směru:

• na pravé straně pole táhne nahoru,
• na levé straně táhne dolů.

V obou případech jde vlastně o to, že se pole mění v prostoru a tyto změny vytvářejí malý vír.

Abychom vír dokázali popsat matematicky, musíme se zamyslet nad tím, co vlastně způsobuje otáčení.

Představme si opět malý čtverec v prostoru. Vektorové pole v něm popisuje například rychlost proudění kapaliny nebo větru. Čtverec je velmi malý, takže pole se v jeho okolí mění jen trochu.

Teď se podívejme na dvě věci.

Nejprve sledujme svislou složku pole $F_y$. Ta říká, jak silně pole tlačí nahoru nebo dolů. Pokud se tato složka mění při pohybu zleva doprava, znamená to, že pravá strana čtverce může být tlačena nahoru jinak než levá. Taková nerovnováha může způsobit otáčení čtverce.

Míru této změny popisuje derivace

$$ \frac{\partial F_y}{\partial x}.$$

Ta říká, jak rychle se svislá složka pole mění při pohybu v horizontálním směru.

Podobný efekt může vzniknout i obráceně. Podívejme se na vodorovnou složku pole $F_x$, která říká, jak silně pole tlačí doleva nebo doprava. Pokud se tato složka mění při pohybu zdola nahoru, může horní strana čtverce tlačit jinak než dolní.

Tuto změnu popisuje derivace

$$\frac{\partial F_x}{\partial y}.$$

A právě kombinace těchto dvou efektů určuje, zda se čtverec začne otáčet.

Rotace pole ve dvou rozměrech je proto dána rozdílem

$$ \frac{\partial F_y}{\partial x} - \frac{\partial F_x}{\partial y}.$$

Tento výraz říká, jak moc pole vytváří lokální otáčení v okolí daného bodu.

Pokud je výsledek nulový, pole v tomto místě žádný vír nevytváří. Pokud je nenulový, pole má vířivý charakter a malé kolečko vložené do pole by se začalo otáčet.

Zatím jsme rotaci definovali pouze pro dvourozměrné pole. Výsledkem byl výraz

$$ \frac{\partial F_y}{\partial x} - \frac{\partial F_x}{\partial y}.$$

Možná jste si všimli jedné zvláštnosti: výsledek je jedno číslo, tedy skalár. To může působit trochu překvapivě, protože rotace přece souvisí s otáčením — a otáčení má vždy nějakou osu.

Důvod je jednoduchý. V čistě dvourozměrném prostoru může vír vzniknout pouze jedním způsobem. Malé kolečko vložené do pole se může otáčet jen kolem osy kolmé k rovině.

Pokud tedy pracujeme v rovině $xy$, je tato osa směrem $z$. Rotace, kterou jsme před chvílí spočítali, je vlastně velikost rotace kolem osy $z$.

Jinými slovy, dvourozměrná rotace je ve skutečnosti jen jedna složka trojrozměrného vektoru rotace.

Jakmile přejdeme do trojrozměrného prostoru, situace se změní. Malé kolečko vložené do pole se může otáčet kolem libovolné osy. Proto už rotaci nemůžeme popsat jedním číslem — potřebujeme vektor, který určí směr i velikost rotace.

Pro vektorové pole

$$ \mathbf{F}(x,y,z) = (F_x, F_y, F_z)$$

je rotace definována jako

$$ \nabla \times \mathbf{F} = \left( \frac{\partial F_z}{\partial y} - \frac{\partial F_y}{\partial z}, \frac{\partial F_x}{\partial z} - \frac{\partial F_z}{\partial x}, \frac{\partial F_y}{\partial x} - \frac{\partial F_x}{\partial y} \right) $$

Každá složka tohoto vektoru popisuje rotaci v jiné rovině prostoru:

• první složka odpovídá rotaci v rovině $yz$,
• druhá složka rotaci v rovině $zx$,
• třetí složka rotaci v rovině $xy$.

A právě tato poslední složka je přesně ten výraz, který jsme odvodili v dvourozměrném případě.

Python ukázka rotace

Abychom si rotaci neukazovali jen na jednom jediném víru, vytvoříme tentokrát složitější vektorové pole, ve kterém budou vedle sebe dva víry s opačným směrem otáčení. Jeden bude mít kladnou rotaci a druhý zápornou.

Díky tomu na jednom obrázku uvidíme nejen proudnice pole, ale také to, že rotace je lokální vlastnost, která se může v různých částech prostoru lišit.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import TwoSlopeNorm

# mřížka
x = np.linspace(-3, 3, 150)
y = np.linspace(-3, 3, 150)
X, Y = np.meshgrid(x, y)

# středy vírů
x1, y1 = -1.5, 0
x2, y2 = 1.5, 0

# vzdálenosti
r1 = (X-x1)**2 + (Y-y1)**2
r2 = (X-x2)**2 + (Y-y2)**2

# víry
U1 = -(Y-y1) * np.exp(-r1)
V1 =  (X-x1) * np.exp(-r1)

U2 =  (Y-y2) * np.exp(-r2)
V2 = -(X-x2) * np.exp(-r2)

# výsledné pole
U = U1 + U2
V = V1 + V2

# rotace
dx = x[1] - x[0]
dy = y[1] - y[0]

dVdx = np.gradient(V, dx, axis=1)
dUdy = np.gradient(U, dy, axis=0)

curl = dVdx - dUdy

# symetrická barevná škála
limit = np.max(np.abs(curl))
norm = TwoSlopeNorm(vmin=-limit, vcenter=0, vmax=limit)

# graf
plt.figure(figsize=(9,7))

plt.contourf(X, Y, curl, levels=80, cmap="RdBu_r", norm=norm)
plt.colorbar(label="Rotation")

plt.streamplot(X, Y, U, V, color="black", density=1.6, linewidth=1)

plt.title("Vector field with two vortices")
plt.xlabel("x")
plt.ylabel("y")
plt.axis("equal")

plt.show()

Vytvoření mřížky bodů

Nejprve si připravíme prostor, ve kterém budeme pole počítat.

x = np.linspace(-3, 3, 150)
y = np.linspace(-3, 3, 150)
X, Y = np.meshgrid(x, y)

Pomocí linspace vytvoříme rovnoměrně rozložené body v intervalu od -3 do 3. Funkce meshgrid z nich sestaví dvourozměrnou síť bodů. Každý bod této sítě představuje jedno místo v rovině, kde budeme počítat směr a velikost vektoru.

Umístění dvou vírů

Potom zvolíme středy dvou vírů:

x1, y1 = -1.5, 0
x2, y2 = 1.5, 0

První vír bude vlevo od středu obrázku, druhý vpravo. Díky tomu budou od sebe dostatečně oddělené a půjde dobře vidět, jak se jejich pole skládají.

Vzdálenost od středu každého víru

Abychom mohli určit, jak silně každý vír působí v různých částech prostoru, spočítáme vzdálenost od jeho středu.

r1 = (X-x1)**2 + (Y-y1)**2
r2 = (X-x2)**2 + (Y-y2)**2

Zde nepoužíváme samotnou vzdálenost, ale její druhou mocninu. To je výhodné, protože se pak dobře kombinuje s exponenciální funkcí.

Konstrukce jednotlivých vírů

Nyní definujeme první a druhý vír zvlášť.

U1 = -(Y-y1) * np.exp(-r1)
V1 =  (X-x1) * np.exp(-r1)

U2 =  (Y-y2) * np.exp(-r2)
V2 = -(X-x2) * np.exp(-r2)

Každý vír má dvě složky:

• U je vodorovná složka pole,
• V je svislá složka pole.

První vír je zkonstruován tak, aby se otáčel jedním směrem, druhý opačným směrem. Exponenciální člen np.exp(-r) zajišťuje, že vír je nejsilnější ve svém středu a směrem ven rychle slábne.

To je důležité, protože jinak by pole působilo stejně silně i ve velmi vzdálených místech a obrázek by nebyl přirozený.

Složení výsledného pole

Celkové pole získáme jako součet obou vírů.

U = U1 + U2
V = V1 + V2

To je velmi důležitá myšlenka, která se v matematice i fyzice objevuje často: složitější pole lze vytvořit jako superpozici jednodušších polí.

Výpočet rotace

Ve dvou rozměrech jsme rotaci definovali jako

$$\frac{\partial F_y}{\partial x} - \frac{\partial F_x}{\partial y}.$$

V našem kódu platí:

• $F_x = U$
• $F_y = V$

Proto počítáme:

dx = x[1] - x[0]
dy = y[1] - y[0]

dVdx = np.gradient(V, dx, axis=1)
dUdy = np.gradient(U, dy, axis=0)

curl = dVdx - dUdy

Funkce np.gradient numericky aproximuje parciální derivace. Jejich rozdílem dostaneme rotaci v každém bodě mřížky.

Výsledkem je skalární pole rotace, které nám říká, kde má pole tendenci roztočit malé kolečko jedním směrem a kde směrem opačným.

Barevná škála se středem v nule

Protože rotace může být kladná i záporná, chceme, aby nula ležela přesně uprostřed barevné škály.

limit = np.max(np.abs(curl))
norm = TwoSlopeNorm(vmin=-limit, vcenter=0, vmax=limit)

To je důležité pro správnou interpretaci obrázku:

• kladná rotace bude mít jednu barvu,
• záporná rotace jinou,
• nulová rotace bude uprostřed.

Vykreslení obrázku

Nejprve vykreslíme barevnou mapu rotace:

plt.contourf(X, Y, curl, levels=80, cmap="RdBu_r", norm=norm)
plt.colorbar(label="Rotation")

Použitá barevná mapa RdBu_r je vhodná pro veličiny, které mohou nabývat kladných i záporných hodnot.

Potom přes ni vykreslíme proudnice vektorového pole:

plt.streamplot(X, Y, U, V, color="black", density=1.6, linewidth=1)

Proudnice ukazují, jak by se v poli pohybovala částice. Díky nim vidíme samotnou strukturu vírů.

Nakonec nastavíme popisky os a poměr stran:

plt.title("Vector field with two vortices")
plt.xlabel("x")
plt.ylabel("y")
plt.axis("equal")

Poměr stran equal je důležitý, aby se víry nedeformovaly.

Jak obrázek interpretovat

Na výsledném obrázku uvidíme dvě důležité vrstvy informace.

První vrstvou jsou proudnice, které ukazují, že se pole v levé a pravé části skutečně otáčí opačnými směry.

Druhou vrstvou je barevná mapa rotace, která ukazuje, že:

• v jedné oblasti je rotace kladná,
• v druhé záporná,
• mezi nimi může být rotace blízká nule.

To je velmi důležité. Rotace není globální vlastnost celého pole, ale lokální vlastnost každého bodu.

Globální vlastnosti pole

Pozorný čtenář si jistě všiml, že jsme zatím mluvili o vlastnostech pole v jednom konkrétním bodě. Gradient, divergence i rotace nám říkají, co se děje v bezprostředním okolí daného místa.

V meteorologii je taková informace velmi užitečná. Díky modelům a měřením můžeme například určit počasí přímo v místě, kde se člověk nachází – jaká je tam teplota, tlak nebo jaký má vítr směr a rychlost – aniž bychom museli tyto hodnoty měřit v každém jednotlivém bodě atmosféry.

Jenže v praxi nás často zajímá ještě něco jiného. Nejen co se děje v jednom bodě, ale kolik „něčeho“ projde určitou oblastí prostoru.

Například:

• kolik vzduchu projde otevřeným oknem,
• kolik elektrického náboje nebo proudu projde vodičem,
• kolik vody proteče potrubím,
• kolik prachových částic zachytí čistička vzduchu.

V těchto případech už nestačí znát hodnotu pole v jednom bodě. Potřebujeme popsat, jak pole prochází určitou plochou nebo objemem. Matematika proto zavádí další důležitý pojem — tok pole.

Tok pole

Doposud jsme se dívali na pole v jednotlivých bodech prostoru. Například jsme mohli zjistit, jaká je rychlost větru nebo jak silné je elektrické pole v konkrétním místě.

V praxi nás ale často nezajímá jen to, co se děje v jednom bodě, ale spíše kolik „něčeho“ projde určitou plochou.

Například:

• kolik vzduchu projde otevřeným oknem,
• kolik vody proteče potrubím,
• kolik elektrického proudu projde vodičem,
• kolik částic zachytí filtr nebo čistička vzduchu.

Ve všech těchto případech nás zajímá tok pole skrz plochu.

Tok pole popisuje, kolik dané veličiny prochází určitou plochou. Pokud máme například pole rychlosti proudění vzduchu, tok nám říká, kolik vzduchu projde danou plochou za jednotku času.

Intuitivně si to můžeme představit jako síť ponořenou do proudu vody. Čím rychleji voda proudí a čím větší je síť, tím více vody skrz ni projde. Naopak pokud síť otočíme šikmo vůči proudu, projde jí méně vody.

Na obrázku vidíme několik ploch (červené desky), které mají stejnou velikost, ale jsou různě natočené vůči směru proudění. Proudění je znázorněno průhlednými válci, které představují proud částic nebo tekutiny procházející prostorem.

Můžeme si představit například proud vzduchu nebo vody. Každý válec znázorňuje malý „sloupec“ proudu, který se pohybuje určitým směrem.

Z obrázku je vidět několik důležitých věcí:

• Spodní vodorovná plocha zachytí největší množství proudu.
• Šikmé plochy zachytí méně proudu.
• Čím více je plocha rovnoběžná se směrem proudění, tím méně proudu skrz ni prochází.

Tok tedy nezávisí jen na velikosti plochy. Záleží také na tom, jak je plocha natočená vůči směru pole.

Musíme proto najít matematický vztah, který tuto závislost popíše.

Jednoduché případy

Začněme nejjednodušší situací.

Pokud pole míří přímo kolmo na plochu, projde skrz plochu maximální množství proudu.

V tomto případě je tok prostě

$$ \Phi = F S $$

kde

• $F$ je velikost pole (například rychlost proudění),
• $S$ je obsah plochy.

Teď si ale představme, že plochu otočíme.

Pokud bude plocha zcela rovnoběžná s prouděním, proud kolem ní jen klouže a skrz plochu nic neprochází. Tok je tedy nulový.

Mezi těmito dvěma extrémy existují všechny možné mezilehlé případy.

Jak popsat natočení plochy

Abychom mohli natočení plochy matematicky popsat, potřebujeme určit směr, který je pro plochu typický.

U každé rovinné plochy můžeme najít směr, který je kolmý na celou plochu. Představ si například šipku, která z plochy „trčí ven“ pod pravým úhlem.

Takové šipce říkáme normála plochy.

Normála tedy není nic složitého. Je to prostě vektor, který je kolmý na danou plochu.

Teď můžeme porovnat dva směry:

• směr pole
• směr normály plochy

Úhel mezi těmito dvěma směry označíme $ \theta $.

Kolmá složka pole

Z fyzikálního hlediska nás nezajímá celá velikost pole, ale jen ta část pole, která skutečně prochází skrz plochu.

Pokud pole míří šikmo, část pole prochází skrz plochu a část pole jen klouže podél ní.

Potřebujeme tedy vzít pouze kolmou složku pole.

Z geometrie víme, že kolmá složka vektoru na daný směr je

$ F \cos\theta $

Proto můžeme tok zapsat jako

$$ \Phi = F S \cos\theta $$

kde

• $F$ je velikost pole,
• $S$ je obsah plochy,
• $\theta$ je úhel mezi polem a normálou plochy.

Elegantní matematický zápis

Výraz

$$ F S \cos\theta $$

není náhodný. V matematice se přesně tento tvar objevuje při určité operaci mezi dvěma vektory.

Této operaci říkáme skalární součin.

Pokud tedy zavedeme vektor plochy $\mathbf{S}$, který

• má směr kolmý k ploše
• a jeho velikost je rovna obsahu plochy

můžeme tok zapsat velmi stručně:

$$\Phi = \mathbf{F} \cdot \mathbf{S}$$

Skalární součin se tedy objevuje přirozeně jako matematický nástroj, který vyjadřuje jak velká část jednoho vektoru míří ve směru druhého.

Co když pole není všude stejné nebo plocha není rovinná?

V ideálním světě je pole všude stejné a plocho je rovina, ale co když chceme popsat reálný svět, kde se velikost pole mění a plocha má od roviny daleko?

Například:

• rychlost proudění vzduchu se v různých místech mění,
• proud vody v potrubí je rychlejší uprostřed než u stěn,
• elektrické pole kolem náboje není všude stejné,
• a plocha, kterou zkoumáme, může být zakřivená (například povrch koule).

V takových situacích už nemůžeme použít jediný jednoduchý vzorec pro celou plochu.

Rozdělení plochy na malé části

Trik, který fyzika a matematika používají, je velmi jednoduchý.

Plochu si rozdělíme na mnoho malých kousků.

Každý z těchto malých kousků je:

• téměř rovinný
• a pole se na něm téměř nemění.

Na tak malém kousku tedy můžeme použít náš předchozí vztah

$$ \Delta \Phi = F \Delta S \cos\theta $$

kde $ \Delta S $ je malý kousek plochy.

Na obrázku vidíme zakřivenou plochu (oranžové těleso). Taková plocha už není rovinná a pole na ní navíc může mít v různých místech různý směr nebo velikost. V takovém případě nemůžeme tok spočítat jedním jednoduchým vzorcem pro celou plochu.

Tok přes celý povrch pak dostaneme tak, že sečteme příspěvky ze všech malých kousků.

Součet všech příspěvků

Pokud bychom plochu rozdělili třeba na tisíce nebo miliony malých kousků, dostali bychom velký součet

$$ \Phi = \Delta \Phi_1 + \Delta \Phi_2 + \Delta \Phi_3 + \cdots $$

Matematika má pro takový součet speciální nástroj.

Když jsou jednotlivé kousky nekonečně malé, součet přechází v integrál.

Tok přes obecnou plochu

Celkový tok pole skrz plochu tedy zapíšeme jako

$$ \Phi = \iint_S \mathbf{F} \cdot d\mathbf{S} $$

Tento zápis znamená:

• rozděl plochu na malé elementy (dS),
• u každého zjisti tok pole,
• a všechny tyto malé příspěvky sečti.

Takový integrál se nazývá plošný integrál vektorového pole.

Formální matematický zápis může na první pohled působit poněkud odstrašujícím dojmem. Ve skutečnosti ale skrývá poměrně jednoduchou myšlenku.

Možná vás napadne, proč se ve vzorci objevují dva znaky integrálu. Důvod je prostý: počítáme tok přes plochu. Plocha je dvojrozměrný objekt, a proto musíme sečíst příspěvky ze dvou směrů zároveň.

Představte si, že plochu rozdělíme na malé čtverečky. Ty můžeme uspořádat do mřížky podobně jako řádky a sloupce v tabulce. Celkový tok pak dostaneme tak, že sečteme příspěvky ze všech těchto malých plošek.

Matematicky se takový součet zapisuje jako dvojný integrál

$$ \iint_S \mathbf{F} \cdot d\mathbf{S} $$

Pokud bychom chtěli tento zápis rozepsat podrobněji, můžeme element plochy vyjádřit pomocí dvou malých změn souřadnic, například $dx$ a $dy$:

$$ \Phi = \iint_S \mathbf{F} \cdot d\mathbf{S} = \int \int \mathbf{F}(x,y) \cdot \mathbf{n} \ dx dy $$

Zápis $dx,dy$ nám připomíná, že vlastně sčítáme malé příspěvky nejprve v jednom směru a potom ve druhém. Dvojný integrál tedy není nic mystického – je to jen kompaktní způsob, jak zapsat součet přes celou plochu.

Jednoduchý integrál sčítá malé úseky čáry. Dvojný integrál sčítá malé plošky. Trojný integrál sčítá malé objemy.

Intuice

Celý postup je vlastně velmi jednoduchý:

1. malý kousek plochy → spočítáme tok
2. všechny malé kousky → sečteme
3. součet nekonečně malých příspěvků → integrál

Python ukázka toku vektorového pole

Předchozí části článku jsme věnovali především intuici. Ukázali jsme si, že tok pole přes plochu lze chápat jako součet malých příspěvků přes malé elementy plochy. Teď si tuto myšlenku vyzkoušíme v praxi pomocí Pythonu.

Budeme zkoumat jednoduché radiální pole, které směřuje od středu prostoru ven. Takové pole se objevuje například ve fyzice:

• u elektrického pole bodového náboje,
• u gravitačního pole bodové hmoty,
• nebo u modelu proudění kapaliny ze zdroje.

Zvolíme kulovou plochu kolem středu pole a spočítáme, kolik pole touto plochou prochází. Program postupuje přesně podle myšlenky, kterou jsme si vysvětlili:

1. povrch koule rozdělíme na mnoho malých plošek,
2. pro každou plošku spočítáme malý příspěvek toku,
3. všechny příspěvky sečteme.

Výsledkem bude numerická aproximace plošného integrálu. Zároveň si pole i plochu zobrazíme, aby bylo jasné, co přesně integrujeme.

Tato ukázka tedy demonstruje dvě věci:

• jak lze tok pole numericky vypočítat,
• a jak souvisí matematický zápis plošného integrálu s jednoduchým algoritmem, který jen sčítá mnoho malých příspěvků.

import numpy as np
import matplotlib.pyplot as plt

# ---------- výpočet toku ----------

R = 1
n_theta = 150
n_phi = 150

theta = np.linspace(0, np.pi, n_theta)
phi = np.linspace(0, 2*np.pi, n_phi)

dtheta = theta[1] - theta[0]
dphi = phi[1] - phi[0]

flux = 0

for t in theta:
    for p in phi:

        x = R*np.sin(t)*np.cos(p)
        y = R*np.sin(t)*np.sin(p)
        z = R*np.cos(t)

        r = np.sqrt(x*x+y*y+z*z)

        # pole
        F = np.array([x,y,z]) / r**3

        # normála koule
        n = np.array([x,y,z]) / r

        # plošný element
        dS = R**2 * np.sin(t) * dtheta * dphi

        flux += np.dot(F, n) * dS

print("Flux:", flux)
print("Expected:", 4*np.pi)


# ---------- vizualizace ----------

fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(projection='3d')

# koule
theta_v = np.linspace(0, np.pi, 40)
phi_v = np.linspace(0, 2*np.pi, 40)

theta_v, phi_v = np.meshgrid(theta_v, phi_v)

x = R*np.sin(theta_v)*np.cos(phi_v)
y = R*np.sin(theta_v)*np.sin(phi_v)
z = R*np.cos(theta_v)

ax.plot_surface(x, y, z, alpha=0.3)

# body pro pole
grid = np.linspace(-1.5,1.5,6)
X,Y,Z = np.meshgrid(grid,grid,grid)

r = np.sqrt(X**2 + Y**2 + Z**2)
r[r==0] = 1

Fx = X/r**3
Fy = Y/r**3
Fz = Z/r**3

ax.quiver(X,Y,Z,Fx,Fy,Fz,length=0.3)

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show()

Rozdělení povrchu koule

Nejprve si připravíme parametrizaci povrchu koule pomocí dvou úhlů.

R = 1

theta = np.linspace(0, np.pi, 150)
phi = np.linspace(0, 2*np.pi, 150)

dtheta = theta[1] - theta[0]
dphi = phi[1] - phi[0]

Povrch koule lze popsat dvěma úhly:

• θ (theta) – úhel od severního pólu k jižnímu
• φ (phi) – úhel kolem svislé osy

Funkce linspace vytvoří rovnoměrně rozložené hodnoty těchto úhlů.
Rozdíly dtheta a dphi pak určují velikost malého elementu plochy.

Tím jsme povrch koule rozdělili na mnoho malých plošek.

Proměnná pro součet toku

Nyní si připravíme proměnnou, do které budeme postupně přičítat malé příspěvky toku.

flux = 0

Každá malá ploška přidá do výsledku malý příspěvek.
Celkový tok pak získáme jejich součtem.

Procházení všech plošek

Teď postupně projdeme všechny kombinace úhlů.

for t in theta:
    for p in phi:

Každá dvojice úhlů odpovídá jednomu malému místu na povrchu koule.

Souřadnice bodu na kouli

Z úhlů spočítáme souřadnice bodu na povrchu koule.

x = R*np.sin(t)*np.cos(p)
y = R*np.sin(t)*np.sin(p)
z = R*np.cos(t)

To je standardní převod ze sférických souřadnic do kartézských.

Každý takový bod představuje střed malé plošky, přes kterou budeme počítat tok.

Vzdálenost od středu

Pro jistotu spočítáme vzdálenost bodu od středu.

r = np.sqrt(x*x + y*y + z*z)

U bodů na kouli je tato vzdálenost rovna poloměru (R), ale výpočet ponecháme obecný.

Definice pole

Nyní definujeme samotné vektorové pole.

F = np.array([x, y, z]) / r**3

Toto pole směřuje radiálně od středu ven.

Podobné pole se objevuje například ve fyzice:

• u elektrického pole bodového náboje
• u gravitačního pole bodové hmoty
• u proudění kapaliny ze zdroje

Velikost pole klesá se vzdáleností.

Normála k povrchu koule

Potřebujeme také určit směr kolmý k povrchu.

n = np.array([x, y, z]) / r

U koule je normála velmi jednoduchá: směřuje od středu ven.

Tento směr je důležitý, protože tok počítáme pouze z části pole, která prochází kolmo skrz plochu.

Velikost plošky

Teď spočítáme obsah malé plošky.

dS = R**2 * np.sin(t) * dtheta * dphi

Tento vzorec vychází z geometrie koule.

Malý element plochy je dán součinem:

• $R^2$
• $\sin(\theta)$
• malých změn úhlů $d\theta$ a $d\phi$

Malý příspěvek toku

Teď spočítáme tok přes jednu malou plošku.

flux += np.dot(F, n) * dS

Operace np.dot(F, n) spočítá skalární součin mezi polem a normálou.

Ten nám říká, jak velká část pole skutečně prochází skrz plochu.

Výsledek pak vynásobíme obsahem malé plošky dS.

Tím získáme malý příspěvek toku.

Celkový tok

Po průchodu všemi ploškami dostaneme celkový tok.

print("Flux:", flux)
print("Expected:", 4*np.pi)

Pro toto pole je známý přesný výsledek:

$\Phi = 4\pi$

Numerický výpočet by měl dát velmi podobné číslo.

To potvrzuje, že součet malých příspěvků skutečně odpovídá plošnému integrálu.

Vizualizace pole a plochy

Ve druhé části programu vykreslíme kouli a vektorové pole.

Koule představuje plochu, přes kterou tok počítáme.
Šipky pak ukazují směr a velikost pole v prostoru.

Díky tomu je vidět, co vlastně integrujeme: pole, které vychází ze středu a prochází povrchem koule.

Jak výsledek interpretovat

Výpočet ukazuje důležitou vlastnost tohoto pole.

Tok přes kouli nezávisí na jejím poloměru.
Pokud bychom kouli zvětšili nebo zmenšili, výsledek by byl stejný.

To znamená, že celkový tok je určen pouze zdrojem pole ve středu.

Právě tato myšlenka stojí za jedním z nejkrásnějších výsledků vektorové analýzy — Gaussovou větou.

Gaussova věta

V předchozí ukázce jsme počítali tok pole přes povrch koule o určitém poloměru. Nabízí se ale několik přirozených otázek.

Co se stane, když změníme poloměr koule? Bude tok větší nebo menší?

A co když místo koule použijeme úplně jiný tvar plochy? Například krychli, válec nebo jiný uzavřený povrch.

Změní se celkový tok pole, nebo zůstane stejný?

Na první pohled by se mohlo zdát, že větší plocha zachytí více pole, a tedy že tok poroste s velikostí nebo tvarem plochy. Ve skutečnosti se ale může stát něco překvapivého.

Pojďme si to vyzkoušet v malém Python experimentu, kde budeme počítat tok pole pro různé plochy a různé velikosti.

import numpy as np

# ------------------------------------------------
# Definice vektorového pole
# ------------------------------------------------
# Radiální pole podobné elektrickému poli náboje

def field(x, y, z):
    r = np.sqrt(x*x + y*y + z*z)
    if r == 0:
        return np.array([0,0,0])
    return np.array([x, y, z]) / r**3


# ------------------------------------------------
# Tok přes kouli
# ------------------------------------------------

def flux_sphere(R, n=200):

    theta = np.linspace(0, np.pi, n)
    phi = np.linspace(0, 2*np.pi, n)

    dtheta = theta[1]-theta[0]
    dphi = phi[1]-phi[0]

    flux = 0

    for t in theta:
        for p in phi:

            x = R*np.sin(t)*np.cos(p)
            y = R*np.sin(t)*np.sin(p)
            z = R*np.cos(t)

            F = field(x,y,z)

            # normála koule
            n_vec = np.array([x,y,z]) / R

            # element plochy
            dS = R**2 * np.sin(t) * dtheta * dphi

            flux += np.dot(F, n_vec) * dS

    return flux


# ------------------------------------------------
# Tok přes krychli
# ------------------------------------------------

def flux_cube(a, n=80):

    flux = 0
    x = np.linspace(-a,a,n)
    y = np.linspace(-a,a,n)
    z = np.linspace(-a,a,n)

    dx = x[1]-x[0]
    dy = y[1]-y[0]
    dz = z[1]-z[0]

    # šest stěn krychle

    for X in x:
        for Y in y:

            # horní stěna
            F = field(X,Y,a)
            flux += np.dot(F,[0,0,1]) * dx*dy

            # dolní stěna
            F = field(X,Y,-a)
            flux += np.dot(F,[0,0,-1]) * dx*dy

    for X in x:
        for Z in z:

            # přední stěna
            F = field(X,a,Z)
            flux += np.dot(F,[0,1,0]) * dx*dz

            # zadní stěna
            F = field(X,-a,Z)
            flux += np.dot(F,[0,-1,0]) * dx*dz

    for Y in y:
        for Z in z:

            # pravá stěna
            F = field(a,Y,Z)
            flux += np.dot(F,[1,0,0]) * dy*dz

            # levá stěna
            F = field(-a,Y,Z)
            flux += np.dot(F,[-1,0,0]) * dy*dz

    return flux


# ------------------------------------------------
# Experiment
# ------------------------------------------------

print("Flux through spheres:")

for R in [0.5, 1, 2, 4]:
    print("R =", R, " -> ", flux_sphere(R))


print("\nFlux through cubes:")

for a in [0.5, 1, 2, 4]:
    print("a =", a, " -> ", flux_cube(a))


print("\nExpected value:", 4*np.pi)

Definice vektorového pole

Nejprve si připravíme samotné pole, jehož tok budeme zkoumat.

import numpy as np

def field(x, y, z):
    r = np.sqrt(x*x + y*y + z*z)
    if r == 0:
        return np.array([0, 0, 0])
    return np.array([x, y, z]) / r**3

Toto pole směřuje v každém bodě radiálně od středu ven. Jeho velikost klesá se vzdáleností, podobně jako u elektrického pole bodového náboje nebo gravitačního pole bodové hmoty.

Zápis

$$ \mathbf{F}(x,y,z) = \frac{(x,y,z)}{r^3} $$

znamená, že směr pole je stejný jako směr vektoru od počátku k bodu ((x,y,z)), ale velikost pole slábne s rostoucí vzdáleností.

U bodu $(0,0,0)$ bychom dělili nulou, proto zde vracíme nulový vektor. V našem experimentu nám to nevadí, protože počátek neleží přímo na integrovacích plochách, ale pouze uvnitř nich.

Výpočet toku přes kouli

Teď vytvoříme funkci, která spočítá tok přes kulovou plochu o zadaném poloměru.

def flux_sphere(R, n=200):

    theta = np.linspace(0, np.pi, n)
    phi = np.linspace(0, 2*np.pi, n)

    dtheta = theta[1] - theta[0]
    dphi = phi[1] - phi[0]

    flux = 0

    for t in theta:
        for p in phi:

            x = R * np.sin(t) * np.cos(p)
            y = R * np.sin(t) * np.sin(p)
            z = R * np.cos(t)

            F = field(x, y, z)

            n_vec = np.array([x, y, z]) / R

            dS = R**2 * np.sin(t) * dtheta * dphi

            flux += np.dot(F, n_vec) * dS

    return flux

V této části děláme přesně to, co jsme si vysvětlili teoreticky.

Povrch koule popíšeme pomocí dvou úhlů:

• $\theta$ určuje polohu od severního pólu k jižnímu,
• $\phi$ určuje otočení kolem svislé osy.

Tím rozdělíme povrch koule na mnoho malých plošek. Pro každou plošku pak:

1. spočítáme její souřadnice,
2. zjistíme hodnotu pole v daném bodě,
3. určíme normálu ke kouli,
4. spočítáme malý příspěvek toku,
5. a ten přičteme k celkovému výsledku.

Zvlášť důležitý je řádek

np.dot(F, n_vec) * dS

Skalární součin np.dot(F, n_vec) určuje, jak velká část pole skutečně prochází skrz plochu. Násobení dS pak zohlední velikost malé plošky.

Celkový tok získáme jako součet všech těchto malých příspěvků.

Výpočet toku přes krychli

Teď uděláme totéž pro jiný tvar plochy, konkrétně pro krychli.

def flux_cube(a, n=80):

    flux = 0
    x = np.linspace(-a, a, n)
    y = np.linspace(-a, a, n)
    z = np.linspace(-a, a, n)

    dx = x[1] - x[0]
    dy = y[1] - y[0]
    dz = z[1] - z[0]

    for X in x:
        for Y in y:

            F = field(X, Y, a)
            flux += np.dot(F, [0, 0, 1]) * dx * dy

            F = field(X, Y, -a)
            flux += np.dot(F, [0, 0, -1]) * dx * dy

    for X in x:
        for Z in z:

            F = field(X, a, Z)
            flux += np.dot(F, [0, 1, 0]) * dx * dz

            F = field(X, -a, Z)
            flux += np.dot(F, [0, -1, 0]) * dx * dz

    for Y in y:
        for Z in z:

            F = field(a, Y, Z)
            flux += np.dot(F, [1, 0, 0]) * dy * dz

            F = field(-a, Y, Z)
            flux += np.dot(F, [-1, 0, 0]) * dy * dz

    return flux

U krychle je situace trochu jiná než u koule, ale základní myšlenka zůstává stejná.

Krychle má šest stěn. Každou z nich rozdělíme na malou mřížku čtverečků. Přes každý takový malý čtvereček spočítáme příspěvek toku a všechny příspěvky sečteme.

Výhodou krychle je, že na každé stěně je normála velmi jednoduchá:

• horní stěna má normálu ((0,0,1)),
• dolní ((0,0,-1)),
• pravá ((1,0,0)),
• levá ((-1,0,0)),
• přední ((0,1,0)),
• zadní ((0,-1,0)).

Proto nemusíme normálu počítat zvlášť pro každý bod. Je na celé stěně stejná.

Spuštění experimentu

Nakonec funkce zavoláme pro různé velikosti koulí a krychlí.

print("Flux through spheres:")
for R in [0.5, 1, 2, 4]:
    print("R =", R, " -> ", flux_sphere(R))

print("\nFlux through cubes:")
for a in [0.5, 1, 2, 4]:
    print("a =", a, " -> ", flux_cube(a))

print("\nExpected value:", 4*np.pi)

Tady testujeme dvě věci:

• co se stane, když změníme velikost koule,
• co se stane, když změníme velikost krychle.

Na konci také vypíšeme očekávanou teoretickou hodnotu

$$ 4\pi \approx 12.566 $$

abychom mohli numerické výsledky porovnat.

Výsledky pro koule

Typický výstup může vypadat například takto:

Vidíme, že při změně poloměru koule se tok prakticky nemění.

To je velmi zajímavý výsledek. Větší koule má sice větší povrch, ale pole je na ní slabší. Menší koule má menší povrch, ale pole je na ní silnější. Tyto dva efekty se navzájem vyruší.

Výsledky pro krychle

Podobně může vypadat výstup pro krychle:

Ani zde se tok prakticky nemění.

To je ještě překvapivější než u koulí. Nezměnili jsme totiž jen velikost plochy, ale i její tvar. Přesto dostáváme téměř stejný výsledek.

Malé rozdíly jsou způsobeny tím, že integrál počítáme numericky, tedy jen přibližně pomocí konečně velké mřížky bodů.

Porovnání s teoretickou hodnotou

Teoretická hodnota je

Numerické výsledky se této hodnotě velmi dobře blíží.

To potvrzuje, že celkový tok tohoto radiálního pole nezávisí na konkrétním tvaru ani velikosti uzavřené plochy, pokud tato plocha obklopuje zdroj v počátku.

Jak výsledek interpretovat

To nejdůležitější není samotné číslo, ale jeho význam.

Experiment ukazuje, že tok není určen tím, zda použijeme kouli, krychli nebo jiný uzavřený povrch. Rozhodující je pouze to, zda plocha uzavírá zdroj pole.

To je přesně myšlenka, která vede ke Gaussově větě.

Můžeme ji zatím formulovat neformálně takto:

Celkový tok uzavřenou plochou závisí na tom, kolik zdrojů pole se nachází uvnitř, a ne na konkrétním tvaru plochy.

Co nám experiment skutečně říká

V experimentu jsme pozorovali zvláštní věc. Tok pole přes uzavřenou plochu se prakticky nezměnil, i když jsme:

• změnili velikost plochy
• změnili její tvar

To naznačuje, že tok nezávisí na geometrii povrchu, ale na něčem jiném. Zkusme se na to podívat intuitivněji. Představme si nádrž s vodou. Uvnitř nádrže je malé čerpadlo, které vodu vytlačuje ven. Kolem čerpadla můžeme postavit různé nádoby:

• kulovou
• krychlovou
• válcovou
• nebo nějaký nepravidelný tvar

Pokud všechny tyto nádoby obklopují čerpadlo, pak množství vody, které za sekundu projde jejich stěnami, bude vždy stejné.

Proč?

Protože voda nikde jinde nevzniká ani nezaniká. Jediný zdroj je čerpadlo. Tok přes stěny nádoby tedy musí být přesně tak velký, kolik vody čerpadlo vytlačí. Teď si představme, že povrch nádoby rozdělíme na mnoho malých plošek. Každou z těchto plošek protéká malé množství vody. Celkový tok pak získáme tak, že tyto malé příspěvky sečteme.

$$ \text{celkový tok} = \text{tok}_1 + \text{tok}_2 + \text{tok}_3 + \dots $$

Přesně tak jsme postupovali i v našem Python programu. Teď si představme, že uvnitř nádoby není jen jedno čerpadlo. Může jich být více. Každé čerpadlo přispívá určitým množstvím vody, které proudí ven. Celkový tok přes stěny nádoby pak bude součet příspěvků všech zdrojů uvnitř.

$$ \text{tok} = \text{zdroj}_1 + \text{zdroj}_2 + \text{zdroj}_3 + \dots $$

V reálném poli ale zdroje nemusí být jen jednotlivé body. Může se stát, že pole vzniká v každém malém kousku prostoru. Představte si například:

• mnoho malých pramínků vody rozptýlených v prostoru
• nebo spoustu malých ventilátorů, které foukají vzduch

Prostor pak můžeme rozdělit na malé objemové buňky. Množství pole vznikající v každé buňce popisuje veličina nazývaná divergence pole. Každá buňka vytváří malý příspěvek pole. Celkový tok přes uzavřenou plochu pak bude součet všech těchto příspěvků:

$$ \text{tok} = q_1 + q_2 + q_3 + \dots $$

kde $q_i$ je „množství pole“, které vzniká v jednotlivých buňkách.

$$ q_i \approx \text{divergence} \times \text{objem}, $$

matematicky by se to vyjádřilo takto:

$$ q_i \approx (\nabla \cdot \mathbf{F})_i \Delta V $$

Tok přes povrch tedy nemusíme počítat jen z pole na povrchu.
Můžeme ho také získat jako součet všech zdrojů pole uvnitř objemu.

Přechod k matematickému zápisu

Pokud bychom prostor rozdělili na velmi malé objemové buňky, mohli bychom psát

$$ \text{tok} = q_1 + q_2 + q_3 + \dots + q_n $$

kde $q_i$ představuje množství pole, které v jednotlivých malých buňkách vzniká nebo zaniká.

Současně můžeme uzavřenou plochu obklopující tento objem rozdělit na mnoho malých plošek. Každou z nich protéká malý příspěvek toku.

Celkový tok přes povrch pak dostaneme jako součet těchto příspěvků:

$$ \text{tok}_1 + \text{tok}_2 + \text{tok}_3 + \dots + \text{tok}_n $$

Musí přitom platit jednoduchá bilance: vše, co uvnitř objemu vznikne nebo zanikne, se musí projevit jako tok přes jeho hranici. Jinými slovy, celkový tok přes uzavřenou plochu je roven součtu všech zdrojů a zániků pole uvnitř objemu.

$$ \text{tok}_1 + \text{tok}_2 + \text{tok}_3 + \dots + \text{tok}_n = q_1 + q_2 + q_3 + \dots + q_n $$

Tato úvaha je velmi obecná. Pole zde může představovat například:

• proudění kapaliny,
• elektrické pole,
• gravitační pole,
• nebo jiné fyzikální pole.

V každé malé objemové buňce přitom množství vznikajícího nebo zanikajícího pole popisuje divergence pole. Divergence tedy udává hustotu zdrojů pole v daném bodě prostoru. Proto se na pravé straně Gaussovy věty objevuje právě divergence – říká totiž, kolik pole v jednotlivých částech objemu vzniká nebo zaniká.

Pokud budeme velikost plošek i objemových buněk zmenšovat stále více, tyto součty přecházejí do matematického zápisu pomocí integrálů. Tím se dostaneme k formální podobě Gaussovy věty.

$$ \iint_S \mathbf{F}\cdot d\mathbf{S} = \iiint_V (\nabla\cdot\mathbf{F}) dV $$

Levá strana představuje tok přes povrch. Pravá strana říká, kolik pole vzniká uvnitř objemu.

Cirkulace pole

Představte si, že jdete pěšky kolem bloku domů za větrného dne.

Vítr fouká určitým směrem. Když jdete po jedné straně bloku, může vám foukat do zad a chůze je příjemná. Na další straně vám fouká přímo proti vám a musíte jít pomaleji. Na jiné straně může foukat téměř kolmo, takže jeho vliv skoro necítíte.

Když obejdete celý blok, můžete si zpětně uvědomit zajímavou věc: někdy vám vítr pomáhal, někdy vás brzdil. Kdybychom dokázali všechny tyto malé příspěvky během cesty sečíst, dostali bychom jedno číslo, které říká, jak silně vítr „obíhá“ kolem této trasy.

Pokud by vítr proudil všude stejným směrem, příspěvky z jednotlivých stran by se navzájem vyrušily a výsledné číslo by bylo téměř nulové.

Ale pokud by se vzduch kolem budov skutečně stáčel do víru, situace by byla jiná. Při obcházení bloku by nám vítr většinu času foukal ve směru naší chůze nebo proti ní, a součet těchto příspěvků by byl výrazně větší.

Abychom tuto myšlenku mohli popsat přesněji, představme si, že naši cestu kolem bloku rozdělíme na mnoho malých úseků.

Na každém z těchto úseků můžeme sledovat, jak moc nám vítr pomáhá nebo škodí. Pokud fouká ve směru naší chůze, jde se nám lépe. Pokud fouká proti nám, musíme vynaložit větší námahu. A pokud fouká kolmo, jeho vliv je malý.

Každý malý úsek cesty tedy přispívá malým „příspěvkem“ k celkovému efektu větru během celé procházky.

Když tyto příspěvky sečteme přes všechny úseky trasy, dostaneme jedno číslo, které popisuje celkový vliv větru při obcházení bloku.

Můžeme si to představit například jako součet práce, kterou musíme během chůze vykonat:

$$ W = w_1 + w_2 + w_3 + \dots + w_n $$

Každý člen $w_i$​ představuje malý příspěvek práce na krátkém úseku cesty.

Velikost tohoto příspěvku závisí především na dvou věcech:

• jak silný je vítr v daném místě,
• jaký je úhel mezi směrem větru a směrem naší chůze.

Pokud vítr fouká přesně ve směru pohybu, jeho vliv je maximální. Pokud fouká proti nám, přispívá záporně. Pokud fouká kolmo, jeho vliv je téměř nulový.

Matematicky tedy každý malý příspěvek závisí na projekci vektoru větru na směr naší chůze.

Když takové příspěvky sečteme podél celé uzavřené trasy, získáme veličinu, která říká, jak silně pole „obíhá“ kolem této křivky.

Zkusme se ale zamyslet nad tím, jak bychom tyto malé příspěvky vlastně počítali.

Představme si, že naši cestu rozdělíme na mnoho velmi krátkých úseků. Každý úsek má určitou délku a směr. Na každém z nich také působí vítr, který má svou velikost a směr.

Na malém úseku cesty nás ale nezajímá celý vítr, ale pouze ta jeho část, která působí ve směru naší chůze. Vítr, který fouká kolmo, nám totiž ani nepomáhá ani nás nebrzdí.

Proto musíme zjistit, jak velká část síly větru leží ve směru našeho pohybu.

Pokud označíme velikost větru $F$, délku malého úseku cesty $r$ a úhel mezi směrem větru a směrem naší chůze jako $\alpha$, pak příspěvek práce na tomto krátkém úseku můžeme zapsat jako

$$ \Delta W = F r \cos(\alpha) $$

Tento vztah říká velmi jednoduchou věc: zajímá nás pouze projekce síly do směru pohybu. Pokud vítr fouká přesně ve směru chůze, $\cos(\alpha)=1$ a jeho vliv je maximální. Pokud fouká proti nám, kosinus je záporný a vítr nás brzdí. Pokud fouká kolmo, kosinus je nulový a vítr na pohyb nemá vliv.

Matematici si všimli, že tento vztah se objevuje velmi často při práci s vektory. Proto pro něj zavedli speciální zápis, kterému říkáme skalární součin.

Předchozí vztah tedy můžeme zapsat stručněji jako

$$ \Delta W = \mathbf{F} \cdot \Delta \mathbf{r} $$

kde

• $\mathbf{F}$ je vektor pole (například vítr),
• $\Delta \mathbf{r}$ je malý úsek naší cesty.

Celkovou práci pak dostaneme jako součet všech těchto malých příspěvků podél celé trasy

$$ W = \mathbf{F}_1 \cdot \Delta \mathbf{r}_1 + \mathbf{F}_2 \cdot \Delta \mathbf{r}_2 + \dots + \mathbf{F}_n \cdot \Delta \mathbf{r}_n $$

Pokud bychom délku těchto úseků zmenšovali stále více, součet jednotlivých příspěvků by se postupně blížil přesnému matematickému zápisu pomocí integrálu.

Tato veličina — která vzniká sčítáním příspěvků pole podél uzavřené křivky — se ve fyzice a matematice nazývá cirkulace pole.

Její formální zápis pak má tvar

$$ \oint_C \mathbf{F}\cdot d\mathbf{r} $$

kde $C$ označuje uzavřenou křivku, podél které se pohybujeme.

Tento integrál tedy vyjadřuje totéž, co jsme si představili intuitivně při obcházení bloku: postupně sčítáme malé příspěvky pole ve směru pohybu po celé trase.

Na první pohled může tento zápis působit trochu tajemně

$$ \oint_C \mathbf{F}\cdot d\mathbf{r} $$

Ve skutečnosti ale neříká nic jiného, než to, co jsme si už představili: sčítáme malé příspěvky pole podél celé uzavřené křivky.

Abychom viděli, jak se takový integrál skutečně počítá, převeďme ho do běžných kartézských souřadnic.

Ve dvourozměrném prostoru můžeme libovolné pole rozložit na jeho složky ve směru os $x$ a $y$

$$ \mathbf{F}(x,y) = (F_x(x,y),,F_y(x,y)) $$

První složka $F_x$ říká, jak silně pole působí ve směru osy $x$, druhá složka $F_y$ ve směru osy $y$.

Když se při obcházení křivky posuneme o malý kousek, změní se naše poloha o malé přírůstky

$$ d\mathbf{r} = (dx, dy) $$

To znamená, že jsme se o trochu posunuli ve směru $x$ a o trochu ve směru $y$.

Malý příspěvek práce jsme si dříve zapsali jako

$$ \mathbf{F}\cdot d\mathbf{r} $$

Pokud ale dosadíme složky obou vektorů, dostaneme obyčejný vztah pro skalární součin

$$ \mathbf{F}\cdot d\mathbf{r} = F_x dx + F_y dy $$

Vidíme tedy, že příspěvek pole se skládá ze dvou částí: jedné způsobené složkou pole ve směru $x$ a druhé způsobené složkou pole ve směru $y$.

Když tyto malé příspěvky postupně sečteme podél celé uzavřené křivky, dostaneme

$$ \oint_C \mathbf{F}\cdot d\mathbf{r} = \oint_C (F_x dx + F_y dy) $$

Tento zápis už není tak tajemný. Říká pouze to, že při obcházení křivky postupně sčítáme dva typy příspěvků: jeden způsobený změnou v $x$ a druhý změnou v $y$.

V praxi se takový integrál často počítá tak, že křivku popíšeme pomocí nějakého parametru. Nejjednodušší je představit si, že po křivce postupně jdeme a náš pohyb popisuje nějaká proměnná $t$. Ta může představovat například čas, ale může to být i jen číslo, které říká, kde se právě na křivce nacházíme.

Poloha bodu na křivce se pak dá popsat dvěma funkcemi

$$ x = x(t), \qquad y = y(t) $$

Ty říkají, jaké souřadnice má bod křivky pro danou hodnotu parametru $t$. Jak se $t$ mění, bod se posouvá po křivce.

Když se $t$ změní o malý kousek $dt$, změní se také naše souřadnice. O malý kousek se posuneme ve směru osy $x$ a o malý kousek ve směru osy $y$. Tyto malé změny můžeme zapsat jako

$$ dx = \frac{dx}{dt} dt, \qquad dy = \frac{dy}{dt} dt $$

To je vlastně jen vyjádření toho, že změna souřadnic závisí na tom, jak rychle se po křivce pohybujeme v jednotlivých směrech.

Dříve jsme si ukázali, že malý příspěvek práce podél křivky je

$$ \mathbf{F}\cdot d\mathbf{r} = F_x dx + F_y dy $$

Teď ale víme, že $dx$ a $dy$ můžeme vyjádřit pomocí parametru $t$. Když je dosadíme, dostaneme

$$ F_x \frac{dx}{dt} dt + F_y \frac{dy}{dt} dt $$

Vidíme tedy, že každý malý příspěvek podél křivky lze vyjádřit jako výraz násobený $dt$. Když tyto příspěvky začneme sčítat, dostaneme integrál podle parametru $t$

$$ \int \left( F_x \frac{dx}{dt} + F_y \frac{dy}{dt} \right) dt $$

Křivkový integrál se tak převede na obyčejný integrál jedné proměnné.

To znamená, že při výpočtu ve skutečnosti děláme přesně to, co jsme si představovali na začátku: postupně procházíme křivku bod po bodu a sčítáme malé příspěvky pole ve směru pohybu.

Symbol

$$ \oint_C \mathbf{F}\cdot d\mathbf{r} $$

tak není nic mystického. Je to jen kompaktní způsob, jak zapsat proces, kdy procházíme celou uzavřenou křivku a sčítáme malé příspěvky pole podél této trasy, přesně tak, jak jsme si to představovali při obcházení bloku ve větru.

Python ukázka cirkulace pole

V předchozí části jsme si ukázali, že cirkulace měří, jak silně pole obíhá kolem uzavřené křivky. Teď si to ukážeme na skutečném fyzikálním příkladu.

Zvolíme magnetické pole kolem přímého vodiče s proudem. Takové pole nevychází ven jako radiální pole bodového náboje, ale naopak vytváří uzavřené kružnice kolem vodiče. To z něj dělá ideální příklad pro výpočet cirkulace.

Budeme tedy počítat, jak velký je součet malých příspěvků pole podél kružnice obepínající vodič. Zároveň celé pole vykreslíme, aby bylo vidět, co vlastně integrujeme.

import numpy as np
import matplotlib.pyplot as plt

# Physical constants
mu0 = 4 * np.pi * 1e-7
I = 1.0

def magnetic_field(x, y):
    r2 = x**2 + y**2
    if r2 == 0:
        return np.array([0.0, 0.0])
    factor = mu0 * I / (2 * np.pi * r2)
    return np.array([-y, x]) * factor

def circulation_circle(R, n=2000):
    theta = np.linspace(0, 2 * np.pi, n, endpoint=False)
    dtheta = 2 * np.pi / n

    circulation = 0.0

    for t in theta:
        x = R * np.cos(t)
        y = R * np.sin(t)

        Bx, By = magnetic_field(x, y)

        dx = -R * np.sin(t) * dtheta
        dy =  R * np.cos(t) * dtheta

        circulation += Bx * dx + By * dy

    return circulation

# Numerical experiment for different radii
radii = [0.5, 1.0, 2.0, 3.0]

print("Numerical circulation around the wire:")
for R in radii:
    value = circulation_circle(R)
    print(f"R = {R:.1f} -> {value:.8e}")

print(f"\nExpected value (mu0 * I): {mu0 * I:.8e}")

# Visualization
fig, ax = plt.subplots(figsize=(8, 8))

grid = np.linspace(-4, 4, 25)
X, Y = np.meshgrid(grid, grid)

U = np.zeros_like(X)
V = np.zeros_like(Y)

for i in range(X.shape[0]):
    for j in range(X.shape[1]):
        bx, by = magnetic_field(X[i, j], Y[i, j])
        U[i, j] = bx
        V[i, j] = by

ax.streamplot(X, Y, U, V, density=1.4, linewidth=1)

# Draw the integration path
R_vis = 2.0
t = np.linspace(0, 2 * np.pi, 400)
xc = R_vis * np.cos(t)
yc = R_vis * np.sin(t)
ax.plot(xc, yc, linestyle="--", linewidth=2, label="Integration path")

# Mark the wire
ax.scatter(0, 0, s=80, label="Wire")

ax.set_title("Magnetic field around a straight current-carrying wire")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_aspect("equal")
ax.legend()
plt.show()

Vysvětlení kódu

Nejprve si připravíme potřebné knihovny a fyzikální konstanty.

import numpy as np
import matplotlib.pyplot as plt

mu0 = 4 * np.pi * 1e-7
I = 1.0

Knihovna numpy nám poslouží pro výpočty s čísly a vektory, matplotlib použijeme pro vykreslení pole.

Proměnná mu0 je permeabilita vakua, tedy fyzikální konstanta, která vystupuje v popisu magnetického pole. Proměnná I představuje elektrický proud ve vodiči. Pro jednoduchost volíme proud roven jednomu ampéru.

Teď nadefinujeme samotné magnetické pole kolem přímého vodiče.

def magnetic_field(x, y):
    r2 = x**2 + y**2
    if r2 == 0:
        return np.array([0.0, 0.0])
    factor = mu0 * I / (2 * np.pi * r2)
    return np.array([-y, x]) * factor

Tato funkce vrací magnetické pole v bodě $(x,y)$.

Nejprve spočítáme

$ r^2 = x^2 + y^2 $

tedy druhou mocninu vzdálenosti od vodiče, který se nachází v počátku soustavy.

Pokud bychom dosadili bod $(0,0)$, dělili bychom nulou. Proto v tomto jediném bodě vracíme nulový vektor. Fyzikálně tam totiž model ideálního nekonečně tenkého vodiče přestává být dobře definovaný.

Výraz

np.array([-y, x])

určuje směr magnetického pole. Tento směr je kolmý na polohový vektor $(x,y)$, takže výsledné pole obíhá kolem vodiče po kružnicích.

Koeficient

mu0 * I / (2 * np.pi * r2)

zajišťuje správnou velikost pole. Ta klesá se vzdáleností od vodiče.

Nyní si vytvoříme funkci, která spočítá součet malých příspěvků pole podél kružnice.

def circulation_circle(R, n=2000):
    theta = np.linspace(0, 2 * np.pi, n, endpoint=False)
    dtheta = 2 * np.pi / n

    circulation = 0.0

Proměnná R určuje poloměr kružnice, po které budeme cirkulaci počítat.

Kružnici popíšeme pomocí úhlu $\theta$, který se mění od $0$ do $2\pi$. Tím rozdělíme celou kružnici na mnoho malých úseků.

Proměnná dtheta udává velikost jednoho malého kroku.

Proměnná circulation bude postupně shromažďovat výsledný součet.

Teď projdeme všechny malé úseky kružnice.

    for t in theta:
        x = R * np.cos(t)
        y = R * np.sin(t)

Každá hodnota t představuje jeden bod na kružnici.

Pomocí známých parametrických rovnic kružnice

$$ x = R\cos t, \qquad y = R\sin t $$

spočítáme souřadnice tohoto bodu.

Tím vlastně říkáme: budeme obcházet kružnici bod po bodu.

V každém bodě zjistíme hodnotu magnetického pole.

        Bx, By = magnetic_field(x, y)

Funkce magnetic_field vrátí dvě složky pole:

• Bx ve směru osy $x$
• By ve směru osy $y$

To je přesně pole, které v daném místě působí.

Teď spočítáme malý posun po kružnici.

        dx = -R * np.sin(t) * dtheta
        dy =  R * np.cos(t) * dtheta

To je velmi důležitý krok.

Nestačí znát jen body kružnice, musíme znát také malý vektor posunu podél křivky. Ten říká, jakým směrem po kružnici právě jdeme.

Protože platí

$$ x = R\cos t, \qquad y = R\sin t $$

jejich derivace podle $t$ jsou

$$ \frac{dx}{dt} = -R\sin t, \qquad \frac{dy}{dt} = R\cos t $$

a po vynásobení malým krokem (d\theta) dostaneme malé změny

$$ dx = -R\sin t d\theta, \qquad dy = R\cos t d\theta $$

Tyto dvě hodnoty dohromady tvoří malý vektor (d\mathbf{r}), tedy malý krok podél křivky.

Teď spočítáme malý příspěvek a přičteme ho k výsledku.

        circulation += Bx * dx + By * dy

Tady se znovu objevuje stejná myšlenka jako dříve u práce a toku.

Nezajímá nás celé pole, ale pouze ta jeho část, která míří ve směru našeho pohybu po křivce.

Proto počítáme výraz

$$ B_x dx + B_y dy $$

což je kartézský zápis skalárního součinu

$$ \mathbf{B}\cdot d\mathbf{r} $$

Pokud pole míří ve směru pohybu, příspěvek je kladný. Pokud míří proti směru, příspěvek je záporný. Pokud je kolmé, příspěvek je nulový.

Součet všech těchto malých příspěvků dává výslednou hodnotu.

Na konci funkce vrátíme spočítanou hodnotu.

    return circulation

Tím dostaneme numerickou aproximaci výrazu

$$ \oint_C \mathbf{B}\cdot d\mathbf{r} $$

tedy součtu malých příspěvků pole podél celé uzavřené kružnice.

Teď provedeme experiment pro několik různých poloměrů.

radii = [0.5, 1.0, 2.0, 3.0]

print("Numerical circulation around the wire:")
for R in radii:
    value = circulation_circle(R)
    print(f"R = {R:.1f} -> {value:.8e}")

print(f"\nExpected value (mu0 * I): {mu0 * I:.8e}")

Zde chceme ověřit velmi zajímavou fyzikální vlastnost.

Spočítáme cirkulaci pro několik různých kružnic kolem stejného vodiče a porovnáme výsledky.

Nakonec vypíšeme také očekávanou teoretickou hodnotu

$$ \mu_0 I $$

kterou předpovídá Ampèrův zákon.

Pokud je vše správně, měly by numerické výsledky vycházet velmi podobně bez ohledu na poloměr kružnice.

Nyní přichází vizualizace samotného pole.

fig, ax = plt.subplots(figsize=(8, 8))

grid = np.linspace(-4, 4, 25)
X, Y = np.meshgrid(grid, grid)

U = np.zeros_like(X)
V = np.zeros_like(Y)

Nejprve vytvoříme obrázek a osy.

Pak si připravíme mřížku bodů v rovině, ve kterých budeme pole počítat. Tato mřížka pokrývá oblast od $-4$ do $4$ v obou směrech.

Pole U a V budou obsahovat složky magnetického pole v jednotlivých bodech mřížky.

Teď pole skutečně spočítáme v každém bodě mřížky.

for i in range(X.shape[0]):
    for j in range(X.shape[1]):
        bx, by = magnetic_field(X[i, j], Y[i, j])
        U[i, j] = bx
        V[i, j] = by

Pro každý bod mřížky zavoláme funkci magnetic_field a uložíme obě složky pole.

Tím dostaneme numerický popis pole v celé zobrazené oblasti.

Samotné pole pak vykreslíme pomocí proudnic.

ax.streamplot(X, Y, U, V, density=1.4, linewidth=1)

Funkce streamplot vykreslí křivky, které sledují směr pole.

V tomto případě by měly vzniknout kružnice kolem středu, protože právě tak vypadá magnetické pole kolem přímého vodiče.

To je velmi důležité: vizualizace okamžitě ukazuje, že pole kolem vodiče skutečně obíhá.

Přes pole ještě vykreslíme kružnici, po které cirkulaci počítáme.

R_vis = 2.0
t = np.linspace(0, 2 * np.pi, 400)
xc = R_vis * np.cos(t)
yc = R_vis * np.sin(t)
ax.plot(xc, yc, linestyle="--", linewidth=2, label="Integration path")

Tato čárkovaná kružnice představuje uzavřenou trasu, po které se sčítají malé příspěvky $\mathbf{B}\cdot d\mathbf{r}$.

Nakonec vyznačíme samotný vodič a nastavíme vzhled obrázku.

ax.scatter(0, 0, s=80, label="Wire")

ax.set_title("Magnetic field around a straight current-carrying wire")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_aspect("equal")
ax.legend()
plt.show()

Bod v počátku představuje vodič, kterým prochází proud.

Poměr stran nastavujeme na equal, aby se kružnice nedeformovaly do elips.

Výsledný obrázek tak ukazuje dvě věci najednou:

• strukturu magnetického pole
• konkrétní uzavřenou křivku, podél které cirkulaci počítáme

Typický výstup může vypadat například takto:

Očekávaná teoretická hodnota je

Vidíme, že cirkulace vychází prakticky stejná pro všechny kružnice.

To je velmi důležitý výsledek. Nezáleží na tom, jak velkou kružnici kolem vodiče zvolíme. Rozhodující je pouze to, že kružnice obepíná vodič s proudem.

Číslo, které jsme spočítali, vyjadřuje celkový „otáčivý efekt“ pole podél celé uzavřené křivky. Při obcházení kružnice jsme v každém malém kroku sečetli příspěvek pole ve směru pohybu. Pokud by pole kolem vodiče neobíhalo, jednotlivé příspěvky by se navzájem vyrušily a výsledná hodnota by byla blízká nule. V našem případě ale pole skutečně vytváří uzavřené kružnice kolem vodiče, takže většina příspěvků má stejný směr a jejich součet dává nenulové číslo. Toto číslo tedy říká, jak silně pole obíhá kolem dané křivky.

Ve fyzice má tento výsledek velmi konkrétní význam: cirkulace magnetického pole kolem uzavřené křivky je přímo úměrná elektrickému proudu, který touto křivkou prochází. Jinými slovy, nenulová cirkulace zde není náhodná vlastnost pole, ale přímý důsledek toho, že ve středu kružnice teče elektrický proud.

Právě tato myšlenka nás přirozeně dovede k dalšímu důležitému vztahu mezi rotací pole uvnitř plochy a cirkulací po jejím okraji, který se nazývá Stokesova věta.

Stokesova věta

Než budeme pokračovat, podívejte se na toto video.

Ve videu jsou krásně vidět experimenty s víry. Nás ale zajímá především jeden konkrétní efekt. Pokud se do víru dostane list, malý kousek papíru nebo třeba pingpongový míček, nezačne se jen unášet proudem jako celek, ale často se začne také otáčet kolem své vlastní osy.

To je velmi důležité pozorování. Ukazuje totiž, že pole v daném místě nedělá jen to, že částici někam posouvá, ale má také schopnost ji lokálně roztáčet.

Právě tuto vlastnost popisuje rotace pole. Pokud je v daném bodě rotace nenulová, malé těleso vložené do proudění má tendenci se otáčet. Pokud je rotace nulová, může být těleso sice unášeno proudem, ale samo se neroztáčí.

Je dobré si uvědomit, že zde mluvíme o lokálním efektu. Nezajímá nás teď celé proudění v nádobě ani celkový tvar víru, ale to, co by se stalo s malým tělesem v jednom konkrétním místě. Rotace tedy říká, jak silně má pole v daném bodě tendenci roztočit malé kolečko, lístek nebo míček.

Teď ale přichází další zajímavá otázka. Pokud bychom nevložili do pole jen jeden malý míček, ale podívali se na celou oblast, například na kruhovou oblast uvnitř víru, jak spolu souvisí lokální otáčení uvnitř této oblasti a celkové obíhání pole po jejím okraji?

Intuitivně bychom čekali, že pokud se uvnitř plochy nachází mnoho malých míst, ve kterých má pole nenulovou rotaci, musí se to nějak projevit i na chování pole podél hranice této oblasti. Jinými slovy: celkový oběh pole po okraji by měl souviset se součtem všech malých lokálních rotací uvnitř.

Pokračujme ještě chvíli v našem myšlenkovém experimentu s vírem.

Představme si, že oblast uvnitř kružnice rozdělíme na mnoho malých plošek, například na malé čtverečky. V každém z nich může mít pole nějakou lokální rotaci. Pokud bychom do každé takové plošky položili malé kolečko nebo drobný lístek, v některých z nich by se začal otáčet rychleji, v jiných pomaleji a někde by se téměř neotáčel.

Každá malá ploška tedy přispívá určitým malým „rotačním efektem“. Pokud bychom chtěli zjistit, jak silně se pole v celé oblasti stáčí, mohli bychom tyto malé příspěvky jednoduše sečíst přes všechny plošky.

Zajímavé ale je, že při takovém rozdělení se většina efektů mezi sousedními buňkami navzájem vyruší. Pohyb pole podél společných hran dvou sousedních plošek se totiž objeví jednou s opačným znaménkem než podruhé. Jinými slovy, co jedna ploška „přidá“, druhá ploška „odečte“.

Pokud bychom tedy sečetli příspěvky ze všech malých buněk uvnitř oblasti, nakonec by zůstaly jen příspěvky podél vnější hranice celé oblasti.

To je velmi překvapivý, ale zároveň hluboce logický výsledek. Celkový efekt všech malých lokálních rotací uvnitř plochy se totiž projeví jako obíhání pole podél jejího okraje.

V našem experimentu s vírem to znamená, že pokud se uvnitř oblasti nachází mnoho malých míst, kde se pole stáčí a roztáčí drobné částice, musí se tento efekt projevit i na tom, jak pole proudí podél hranice oblasti. Jinými slovy, lokální otáčení uvnitř a celkové obíhání po okraji jsou dvě různé stránky stejného jevu.

Pokud bychom tuto úvahu zapsali pro konečný počet malých buněk, dostali bychom vztah ve tvaru součtu. Představme si, že celou plochu rozdělíme na malé plošky $S_1, S_2, S_3, \dots , S_n$. V každé z nich může mít pole určitou lokální rotaci.

Každá ploška pak přispívá malým příspěvkem, který je úměrný rotaci pole v daném místě a velikosti plošky. Pro konečný počet buněk bychom tedy mohli psát přibližně

$$ \text{obíhání po okraji} \approx (\text{rotace}_1 \Delta S_1) + (\text{rotace}_2 \Delta S_2) + \dots + (\text{rotace}_n \Delta S_n). $$

Zároveň ale víme, že obíhání po okraji jsme už dříve počítali jako součet malých příspěvků pole podél uzavřené křivky

$$ \mathbf{F}_1 \cdot \Delta \mathbf{r}_1 + \mathbf{F}_2 \cdot \Delta \mathbf{r}_2 + \dots + \mathbf{F}_m \cdot \Delta \mathbf{r}_m.$$

Vidíme tedy, že máme dva různé způsoby, jak popsat stejný jev:

• součet malých příspěvků podél okraje oblasti,
• součet malých rotačních efektů uvnitř celé plochy.

Pokud tedy budeme popisovat stejnou situaci dvěma různými způsoby, můžeme napsat vztah mezi těmito součty

$$ \mathbf{F}_1 \cdot \Delta \mathbf{r}_1 + \mathbf{F}_2 \cdot \Delta \mathbf{r}_2 + \dots + \mathbf{F}_m \cdot \Delta \mathbf{r}_m (\text{rotace}_1 \Delta S_1) + (\text{rotace}_2 \Delta S_2) + \dots + (\text{rotace}_n \Delta S_n). $$

Levá strana představuje součet malých příspěvků pole podél uzavřené křivky, zatímco pravá strana představuje součet malých rotačních efektů uvnitř celé plochy, kterou tato křivka ohraničuje.

Jak budeme buňky zmenšovat a jejich počet zvětšovat, tyto součty se budou stále více zpřesňovat. V limitě, kdy jsou plošky i úseky křivky nekonečně malé, přecházejí součty v integrály.

Součet po okraji se změní v křivkový integrál

$$ \oint_C \mathbf{F}\cdot d\mathbf{r}, $$

zatímco součet rotačních příspěvků přes plochu přejde v plošný integrál

$$ \iint_S (\nabla \times \mathbf{F}) \cdot \mathbf{n}, dS. $$

Teprve v tomto okamžiku můžeme celý vztah zapsat kompaktně jako

$$ \oint_C \mathbf{F}\cdot d\mathbf{r} = \iint_S (\nabla \times \mathbf{F}) \cdot \mathbf{n} dS.$$

Tento vztah říká, že celkové obíhání pole podél okraje oblasti je dáno součtem všech malých rotací pole uvnitř této oblasti. V matematice a fyzice se tento vztah nazývá Stokesova věta.

Python experiment a ověření Stokesovy věty na speciálním případě

Teď jsme si vlastně zavedli určitou větu – myšlenku, že obíhání pole po okraji oblasti souvisí se součtem rotací uvnitř této oblasti. Abychom získali důvěru v to, že tato úvaha dává smysl, ukážeme si ji na konkrétním příkladu.

Je fér říct, že historicky to takto neprobíhalo. Já jsem zde udělal malou mentální zkratku. V matematice a fyzice se obecné věty obvykle nerodí tak, že je někdo napíše rovnou v konečné podobě. Mnohem častěji lidé nejprve počítají velké množství konkrétních případů. Postupně si začnou všímat, že se v různých situacích opakují podobné vzorce a vztahy. Teprve poté někdo přijde s myšlenkou, že by za těmito jednotlivými případy mohl stát obecnější princip.

Tento postup je mimochodem velmi blízký tomu, jak funguje náš mozek. Po milionech let evoluce jsme velmi dobří v rozpoznávání vzorců. V minulosti nám to pomáhalo přežít v přírodě – například rozpoznat stopy zvířat nebo změny počasí. Dnes stejný mechanismus používáme při řešení složitých vědeckých a technických problémů.

Uděláme tedy něco podobného. Místo okamžitého přijetí obecné věty si nejprve vezmeme jednoduché vířivé pole a numericky na něm ověříme, že součet příspěvků po okraji oblasti opravdu odpovídá součtu lokálních rotací uvnitř plochy. Na tomto speciálním případu tak uvidíme, že naše úvaha skutečně funguje.

import numpy as np
import matplotlib.pyplot as plt

def field(x, y):
    return np.array([-y, x])

def curl(x, y):
    return 2.0


def circulation_on_circle(R=1.0, n=2000):
    theta = np.linspace(0, 2 * np.pi, n, endpoint=False)
    dtheta = 2 * np.pi / n

    total = 0.0

    for t in theta:
        x = R * np.cos(t)
        y = R * np.sin(t)

        Fx, Fy = field(x, y)

        dx = -R * np.sin(t) * dtheta
        dy =  R * np.cos(t) * dtheta

        total += Fx * dx + Fy * dy

    return total


def curl_flux_through_disk(R=1.0, n=400):
    x = np.linspace(-R, R, n)
    y = np.linspace(-R, R, n)
    dx = x[1] - x[0]
    dy = y[1] - y[0]

    total = 0.0

    for X in x:
        for Y in y:
            if X**2 + Y**2 <= R**2:
                total += curl(X, Y) * dx * dy

    return total


R = 1.5

circulation_value = circulation_on_circle(R=R)
curl_flux_value = curl_flux_through_disk(R=R)
expected_value = 2 * np.pi * R**2

print(f"Circle radius: {R}")
print(f"Circulation along boundary: {circulation_value:.6f}")
print(f"Sum of curl over disk:      {curl_flux_value:.6f}")
print(f"Expected value:             {expected_value:.6f}")


fig, ax = plt.subplots(figsize=(8, 8))

grid = np.linspace(-2.2, 2.2, 25)
X, Y = np.meshgrid(grid, grid)

U = -Y
V = X

ax.streamplot(X, Y, U, V, density=1.4, linewidth=1)

t = np.linspace(0, 2 * np.pi, 400)
xc = R * np.cos(t)
yc = R * np.sin(t)
ax.plot(xc, yc, "--", linewidth=2, label="Boundary curve")

disk = plt.Circle((0, 0), R, fill=False, linewidth=2, alpha=0.5)
ax.add_patch(disk)

ax.set_title("Vector field and boundary of the integration region")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_aspect("equal")
ax.legend()
plt.show()

Nejprve si nadefinujeme samotné pole a jeho rotaci.

def field(x, y):
    return np.array([-y, x])

def curl(x, y):
    return 2.0

Pole $( -y, x )$ je jednoduché vířivé pole. V každém bodě směřuje tečně kolem středu, takže jeho proudnice tvoří kružnice. Je to velmi vhodný model, protože je dostatečně jednoduchý na výpočet a zároveň na něm dobře vidíme vztah mezi lokální rotací uvnitř oblasti a obíháním po jejím okraji.

Funkce curl vrací hodnotu rotace pole. V tomto případě je rotace v každém bodě stejná a rovná dvěma. To znamená, že pole má v celé oblasti stejnou lokální tendenci roztáčet malé kolečko.

Potom vytvoříme funkci, která spočítá součet malých příspěvků pole podél kružnice.

def circulation_on_circle(R=1.0, n=2000):
    theta = np.linspace(0, 2 * np.pi, n, endpoint=False)
    dtheta = 2 * np.pi / n

    total = 0.0

Proměnná R určuje poloměr kružnice, po které budeme obcházet okraj oblasti. Proměnná n určuje, na kolik malých úseků tuto kružnici rozdělíme. Čím větší n, tím přesnější bude výsledek.

Pole theta obsahuje jednotlivé úhly, které odpovídají bodům na kružnici. Proměnná dtheta udává velikost jednoho malého kroku. Proměnná total bude postupně shromažďovat celkový součet příspěvků.

Teď budeme kružnici obcházet bod po bodu.

    for t in theta:
        x = R * np.cos(t)
        y = R * np.sin(t)

        Fx, Fy = field(x, y)

Pomocí parametrických rovnic kružnice spočítáme souřadnice aktuálního bodu. V každém z nich pak zjistíme hodnotu pole. Tím získáme vektor, který v daném místě popisuje směr a velikost pole.

Pak potřebujeme určit malý posun po křivce.

        dx = -R * np.sin(t) * dtheta
        dy =  R * np.cos(t) * dtheta

Tento krok je velmi důležitý. Nestačí znát jen polohu na kružnici, musíme znát také směr, kterým se po ní právě pohybujeme. Hodnoty dx a dy představují malý krok podél křivky. Jsou získané z derivací parametrických rovnic kružnice.

Právě tento malý vektor pohybu hraje stejnou roli jako $d\mathbf{r}$ v křivkovém integrálu.

Nyní spočítáme malý příspěvek a přičteme ho k celkovému součtu.

        total += Fx * dx + Fy * dy

Tady se znovu objevuje stejná myšlenka jako u práce a cirkulace. Nezajímá nás celé pole, ale pouze ta jeho část, která míří ve směru našeho pohybu po křivce. Proto počítáme výraz

$$ F_x dx + F_y dy,$$

což je kartézský zápis skalárního součinu

$$ \mathbf{F}\cdot d\mathbf{r}. $$

Pokud pole míří ve směru pohybu, příspěvek je kladný. Pokud míří proti směru, je záporný. Pokud je kolmé, nepřispěje téměř vůbec.

Na konci funkce vrátíme výslednou hodnotu.

    return total

Tím dostaneme numerickou aproximaci cirkulace pole podél kružnice.

Druhá část kódu počítá součet lokálních rotací uvnitř celé plochy.

def curl_flux_through_disk(R=1.0, n=400):
    x = np.linspace(-R, R, n)
    y = np.linspace(-R, R, n)
    dx = x[1] - x[0]
    dy = y[1] - y[0]

    total = 0.0

Tentokrát už neprocházíme okraj, ale celou oblast uvnitř kruhu. Nejprve si vytvoříme mřížku bodů v čtverci od (-R) do (R). Proměnné dx a dy udávají velikost jedné malé buňky této mřížky. Proměnná total bude opět shromažďovat součet všech příspěvků.

Pak projdeme všechny body mřížky.

    for X in x:
        for Y in y:
            if X**2 + Y**2 <= R**2:
                total += curl(X, Y) * dx * dy

Nejprve testujeme, zda daný bod leží uvnitř kruhu. Pokud ano, přičteme příspěvek dané malé plošky. Ten je dán jako hodnota rotace v daném místě krát obsah malé buňky $dx,dy$.

Tím vlastně děláme přesně to, co jsme si před chvílí popsali intuitivně: bereme malé rotační efekty uvnitř plochy a sčítáme je přes celou oblast.

Na konci opět vrátíme výsledek.

    return total

Tím získáme numerickou aproximaci plošného integrálu rotace přes celý kruh.

Nyní oba výpočty porovnáme.

R = 1.5

circulation_value = circulation_on_circle(R=R)
curl_flux_value = curl_flux_through_disk(R=R)
expected_value = 2 * np.pi * R**2

print(f"Circle radius: {R}")
print(f"Circulation along boundary: {circulation_value:.6f}")
print(f"Sum of curl over disk:      {curl_flux_value:.6f}")
print(f"Expected value:             {expected_value:.6f}")

Zvolíme poloměr kružnice a spočítáme:

• součet příspěvků po okraji,
• součet rotací uvnitř plochy,
• a také teoretickou hodnotu.

Protože rotace je v našem případě konstantní a rovná dvěma, celkový součet přes kruh musí být

$$ 2 \cdot \text{obsah kruhu} = 2\pi R^2. $$

To je třetí způsob, jak si můžeme správnost výsledku ověřit.

Nakonec pole vykreslíme.

fig, ax = plt.subplots(figsize=(8, 8))

grid = np.linspace(-2.2, 2.2, 25)
X, Y = np.meshgrid(grid, grid)

U = -Y
V = X

Nejprve si připravíme obrázek a mřížku bodů, ve kterých budeme pole zobrazovat. Hodnoty U a V představují složky vektorového pole v jednotlivých bodech.

Samotné pole pak vykreslíme pomocí proudnic.

ax.streamplot(X, Y, U, V, density=1.4, linewidth=1)

Proudnice krásně ukážou, že pole skutečně obíhá kolem středu. To je důležité, protože právě takové chování očekáváme od vířivého pole.

Přes obrázek ještě vykreslíme hranici oblasti.

t = np.linspace(0, 2 * np.pi, 400)
xc = R * np.cos(t)
yc = R * np.sin(t)
ax.plot(xc, yc, "--", linewidth=2, label="Boundary curve")

Tato čárkovaná kružnice představuje uzavřenou křivku, po které jsme počítali součet příspěvků pole.

Nakonec zvýrazníme i samotnou oblast a nastavíme vzhled obrázku.

disk = plt.Circle((0, 0), R, fill=False, linewidth=2, alpha=0.5)
ax.add_patch(disk)

ax.set_title("Vector field and boundary of the integration region")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_aspect("equal")
ax.legend()
plt.show()

Kružnice tak na obrázku hraje dvojí roli. Je to jednak okraj, po kterém počítáme cirkulaci, a zároveň ohraničuje plochu, přes kterou sčítáme rotaci pole.

To je na celé ukázce to nejdůležitější. Jeden a tentýž jev zde počítáme dvěma různými způsoby:

• jednou jako součet malých příspěvků po hranici,
• podruhé jako součet malých rotačních efektů uvnitř celé plochy.

Proč je to vlastně tohle všechno důležité

Tohle by mělo navázat jak na článek, tak na tvůj seriál. Navrhuji text v tomto duchu:

Možná vás při čtení napadla otázka: proč jsme vlastně prošli celou tuto teorii? Gradient, divergence, rotace, tok pole, cirkulace, Gaussova a Stokesova věta mohou na první pohled působit jako čistě matematické konstrukce.

Ve skutečnosti ale tvoří základní jazyk pro popis fyzikální reality.

Velká část světa kolem nás se totiž popisuje pomocí polí. Teplota v atmosféře, rychlost proudění vzduchu, tlak v kapalině, elektrické a magnetické pole, koncentrace chemických látek nebo například hustota energie v prostoru. Ve všech těchto případech má každé místo v prostoru svou hodnotu a tyto hodnoty se navzájem ovlivňují.

Abychom takové systémy dokázali pochopit a počítat, potřebujeme nástroje, které popisují jak se pole mění, kde vzniká, kde zaniká a kde se stáčí do vírů. Přesně k tomu slouží pojmy, které jsme si v tomto článku postupně zavedli.

Tyto myšlenky se neobjevují jen v učebnicích matematiky. Jsou přímo zabudované v rovnicích, které popisují svět kolem nás. Najdeme je například v elektromagnetismu, mechanice tekutin, meteorologii, akustice nebo v mnoha oblastech chemického a procesního inženýrství.

Právě díky nim dnes dokážeme vytvářet numerické modely reality. Počítače simulují proudění vzduchu kolem letadel, šíření elektromagnetických vln v anténách, proudění vody v oceánech nebo transport látek v chemických reaktorech. V jádru těchto simulací se stále objevují stejné principy: prostor se rozdělí na velké množství malých buněk a v každé z nich se počítají malé příspěvky polí. Součty těchto příspěvků pak aproximují integrály, které popisují fyzikální zákony.

Právě zde vstupuje do hry Python.

Tento článek je součástí série o modelování reality pomocí Pythonu. Na začátku série jsem napsal větu:

„Zapomeňte na e-shopy. Python umí modelovat realitu.“

A přesně o tom to je. Python dnes není jen jazyk pro webové aplikace nebo práci s daty. Je to nástroj, který umožňuje simulovat fyzikální procesy, modelovat složité systémy a experimentovat s matematickými modely světa kolem nás.

V dalších článcích této série se proto postupně posuneme od teorie k praxi a ukážeme si reálné modely a simulace, kde se všechny tyto myšlenky skutečně používají.

Na konci článku se můžeme na chvíli zastavit a podívat se na celý obraz z nadhledu. Během čtení jsme postupně zavedli několik pojmů, které na první pohled mohou působit nesourodě. Ve skutečnosti ale tvoří velmi kompaktní systém pro popis polí v prostoru.

Gradient říká, jak rychle se skalární pole mění a kterým směrem roste nejrychleji.

Divergence popisuje, zda pole v daném místě vzniká nebo zaniká – tedy zda se z bodu „rozbíhá“ nebo naopak „stahuje“.

Rotace (curl) ukazuje, zda má pole v okolí bodu tendenci stáčet se do víru.

Tyto lokální vlastnosti pole můžeme následně sčítat přes větší oblasti prostoru.

Tok pole popisuje, kolik pole projde danou plochou.

Cirkulace popisuje, jak silně pole „obíhá“ podél uzavřené křivky.

Nakonec přicházejí dvě velmi hluboké věty, které tyto lokální a globální popisy propojují.

Gaussova věta říká, že tok pole přes uzavřený povrch souvisí s divergencí pole uvnitř objemu.

Stokesova věta říká, že cirkulace pole podél uzavřené křivky souvisí s rotací pole uvnitř plochy.

Obě věty vyjadřují stejnou základní myšlenku: globální chování pole na hranici oblasti je určeno tím, co se děje uvnitř této oblasti.

Právě tato myšlenka stojí v pozadí velké části moderní fyziky i numerických simulací. A také proto jsme si ji v tomto článku krok za krokem vybudovali — od intuice a jednoduchých součtů až po integrální zápisy.