Matematický aparát pro modelování fyzikálních polí v příkladech

Technická univerzita v Liberci, Fakulta mechatroniky, informatiky a mezioborových studií

Lukáš Zedek a kolektiv

Matematický aparát pro modelování fyzikálních polí v příkladech

L iberec, 2015

recenzent: doc. RNDr. Edita Kolářová, Ph.D., © Lukáš Zedek, Miloslav Košek, Milan Hokr, Petr Rálek ISBN 978-80-7494-268-6

Za podíl na přípravě vydání této knihy děkuji celému autorskému kolektivu stejně tak jako lidem, kteří mi pomohli svou kontrolou obsahové a jazykové správnosti textu. Jmenovitě bych rád poděkoval: Miloslavu Koškovi, Milanu Hokrovi, Petru Rálkovi, Janu Šemberovi, a Martinu Plešingerovi. Za poskytnutí možnosti připravit tuto publikaci si zaslouží speciální dík prof. Aleš Richter.

Lukáš Zedek

3

Obsah Značení fyzikálních veličin

7

Značení matematických operací a symbolů

8

Seznam zkratek

9

1 Úvod 2 Fyzikální veličiny a zákony 2.1 Souřadná soustava . . . . 2.2 Fyzikální veličiny . . . . . 2.3 Matematický popis . . . . 2.4 Materiálové prostředí . . .

11 . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

13 13 14 15 16

3 Vektorová algebra 19 3.1 Vektorový prostor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2 Vektory ve fyzice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.3 Některé užitečné vztahy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4 Vektorová analýza 4.1 SageMath: snáze, rychleji, přesněji . 4.2 Operátor nabla . . . . . . . . . . . . 4.3 Skalární pole . . . . . . . . . . . . . 4.4 Vektorové pole . . . . . . . . . . . . 4.5 Důležité vztahy vektorové analýzy . 4.6 Matematický popis vektorového pole

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

27 27 27 27 29 49 52

5 Základy tenzorového počtu 5.1 Tenzor elektrické vodivosti . . . . . . 5.2 Tenzory v mechanice pružných těles 5.3 Vektory radiální a axiální . . . . . . 5.4 Tenzory vyššího řádu . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

55 55 57 64 65

6 Analogie 67 6.1 Elektromechanická analogie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 6.2 Analogie mezi transportními jevy . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 6.3 Poissonova a Laplaceova rovnice . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 7 Postupy řešení matematických modelů 7.1 Klasifikace rovnic a jejich použití . . . . 7.2 Analytické řešení rovnic a soustav rovnic 7.3 Numerické řešení . . . . . . . . . . . . . 7.4 Metody zpracování naměřených dat . . . 8 Závěr

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

73 . 73 . 80 . 88 . 105 119

5

Summary The book describes mathematical tools for the description and modelling as well as prediction of an evolution of quantities in natural and technological systems. Necessary theoretical knowledge basis is illustrated by numerous examples with practical application. The reader will get a chance to practice and brush up well-established analytical and numerical methods for a solution of mathematical models. To make the solution as simple as possible, we introduce selected opensource (free) solvers and program tools.

6

Značení fyzikálních veličin − → B c C − → D − → E ϵ ϵ0 f − → F φ φ φ(x, t) φ(x, y, z) G − → G γ i(t), I Iˆ0 − → i k Kv m ω ω ωij p Pˆ − → p Q, q, q(t) − → qv r R ρ ρ − → → − r,− r A, → rB s S Sij t T Tij u(t), U ˆ0 U Uo , U v − → u − → v W, A ZˆR , ZˆL , ZˆC Zˆ

vektor magnetické indukce rychlost šíření vlny, resp. rychlost světla kapacita kondenzátoru vektor elektrické indukce intenzita elektrického pole absolutní permitivita permitivita vakua frekvence kmitání/vlnění vektor síly fáze kmitání/vlnění elektrický potenciál piezometrická výška (potenciál) skalární potenciál elektrická vodivost součástky intenzita potenciálového pole tenzor (měrné) elektrické vodivosti okamžitá hodnota elektrického proudu ve vodiči, elektrický proud fázor elektrického proudu prostorová proudová hustota tuhost pružiny materiálová konstanta definující transportní“ vlastnosti prostředí ” hmotnost úhlová rychlost tenzor pootočení složka tenzoru pootočení tlak komplexní výkon hybnost elektrický náboj matematicky definovaný tok“ veličiny v ” mechanický odpor elektrický odpor objemová hustota volného náboje hustota látky polohový vektor dráha tenzor deformace složka tenzoru deformace čas tenzor mechanického napětí složka tenzoru mechanického napětí elektrické napětí fázor napětí amplituda napětí vektor posunutí rychlost, pole rychlostí mechanická práce impedance rezistoru, induktoru a kapacitoru celková impedance 7

Značení matematických operací a symbolů − → 0 At An A x , Ay , Az − →→ A(− r) det(D) det(D)1,j − → div A − → → d− r,dl − → dS dV dα dΦ δij ∆ ∂Ai ∂xj E (N )

− →, − → − → e 1 e2 , e3 ∃ f ′ (. . .), f (k) (. . .) ϕ(s) Φ Φn grad φ(x, y, z) Γ∫ ∫ ∫ , ,

nulový prvek vektorového prostoru − → tečná složka vektoru A − → normálová složka vektoru A − → složky vektoru A vektorové pole determinant matice D subdeterminant matice D − → divergence vektoru A element dráhy, element křivky element orientované plochy element objemu V dílčí úhel dílčí tok vektoru plochou Kroneckerovo delta Laplaceův operátor parciální derivace složky vektoru N-rozměrný Eukleidovský prostor jednotkové směrové vektory souřadných os kvantifikátor - existuje prvek . . .“ ” derivace funkce f (. . .) Laplaceův obraz funkce φ(t) tok vektoru plochou normovaný tok vektoru plochou gradient skalární funkce hranice oblasti Ω v prostoru integrál po křivce, na ploše, v objemu

(K) (S) (V )

j ∀ λi λi,j L [φ(x, t)] L−1 [Φ0 ] − → n Ω − → ∇ − → ∇φ(x, y, z) − → ∇φ(x1 , x2 , x3 ) → − → − ∇·A − − → → ∇×A ρ(T) → − rot A − → t− → u = u − →·− → v) (u− → → u,− v − → → − u × v] [− → → − u, v

komplexní jednotka kvantifikátor - pro každý prvek . . .“ ” i-té vlastní číslo matice řešení charakteristické rovnice Laplaceova transformace funkce φ(x, t) inverzní Laplaceova transformace obrazu Φ(s) normálový vektor podoblast prostoru vektorový operátor nabla gradient skalární funkce gradient skalární funkce − → divergence vektoru A − → rotace vektoru A spektrální poloměr matice T − → rotace vektoru A tečný vektor norma vektoru skalární součin skalární součin vektorový součin vektorový součin 8

Značení matematických operací a symbolů X, Y, Z X1 , X 2 , X 3 x 1 , x2 , x 3 Vij

označení souřadných os označení souřadných os označení souřadných os prvek na pozici v i-tém řádku a j-tém sloupci matice V

Seznam zkratek AE (AR) DAE (DAR) LS FEM (MKP) ODE (ODR) PCA PDE (PDR) SLAR SVD TLS

Algebraic Equation (Algebraická Rovnice) Differential Algebraic Equation (Diferenciálně Algebraická Rovnice) Least Squares (Metoda Nejmenších Čtverců) Finite Element Method (Metoda Konečných Prvků) Ordinary Differential Equation (Obyčejná Diferenciální Rovnice) Principal Component Analysis (Analýza Hlavních Komponent) Partial Differential Equation (Parciální Diferenciální Rovnice) Soustava Lineárních Algebraických Rovnic Singular Value Decomposition (Singulární Rozklad Matice) Total Least Squares (Úplná Metoda Nejmenších Čtverců)

9

Předmluva Teoretické fyzikální a technické publikace využívají relativně složitý specifický matematický aparát, který je pro běžného čtenáře zpočátku málo srozumitelný. To se týká zejména teorie elektromagnetického pole, viz např. kvalitní českou učebnici [3] nebo rozsáhlou příručku [18]. Tímto oborem se zabývá i řada skript, např. skriptum [11]. Lze odkázat i na rozsáhlou cizojazyčnou literaturu, jako reprezentanta uvádíme prakticky zaměřenou publikaci [5]. Výklad teorie elektromagnetického pole je nemyslitelný bez rozsáhlého využití vektorové analýzy, s níž bývají začátečníci seznámeni jen okrajově. Autoři fyzikálně či technicky zaměřených monografií u čtenáře předpokládají dobrou znalost použitého matematického aparátu a oprávněně se zaměřují na vysvětlení fyzikálního významu, či technického využití výsledků tímto aparátem získaných. Pro úplnost je nutno dodat, že např. v knize [3] je v příloze uveden přehled používaných základních vztahů vektorové algebry a analýzy. Existují ovšem i takové monografie, které ilustrují využití matematického aparátu na fyzikálních příkladech. Pro bližší seznámení se s širokým spektrem matematických nástrojů je čtenáři možno doporučit knihu [6] nebo rozsáhlé přehledy užité matematiky [17] od prof. Rektoryse. Dnešní přístup k výzkumu však vyžaduje okamžité výsledky a tudíž nepodporuje důkladné teoretické studium, které zabere mnoho času a jehož výsledky se výrazně projeví až po letech. Důkladné studium knih o matematickém aparátu se zdánlivě jeví jako ztráta času. Kromě toho čtenář obvykle potřebuje znát podrobně jen tu část z těchto monografií, která je využitelná pro jeho aplikovaný výzkum. Proto jsme připravili stručnou příručku, která se soustřeďuje jen na vybrané části matematického aparátu, které se jeví jako nezbytné pro naše nové, mladé kolegy na fakultě mechatroniky Technické univerzity v Liberci. Příručka se snaží jednoduše, srozumitelně a relativně přesně popsat speciální matematické metody používané ve fyzice, vzápětí upozorňuje na jejich fyzikální význam a na možné technické využití. V publikaci tedy není zdaleka jen čistá matematika, která by od čtení příručky mohla leckoho odradit. Naopak, kniha je na mnoha místech doplněna fyzikou, tedy fyzikálním významem a aplikacemi odvozených nebo vysvětlených vztahů. Ke zvýšení využitelnosti publikace by měly dopomoci jednoduché ilustrativní příklady, kterými jsou proloženy jednotlivé kapitoly. Na základě dlouholetých pedagogických zkušeností jsme do příručky zařadili i části, které s matematickým aparátem přímo nesouvisí, ale týkají se témat ve kterých mnozí lidé často chybují. Jedná se např. o druhou kapitolu, Fyzikální veličiny a zákony. Dále je zařazena kapitola o tenzorech, v níž se důkladněji probírají tenzory druhého řádu se zaměřením správně pochopit význam tenzoru elastického napětí a deformace. To považujeme za důležité, poněvadž elasticita se v mechatronice stále více využívá. V závěrečných kapitolách jsou nastíněny možné aplikace matematického aparátu. V kapitole 6 je zmínka o využití analogií mezi rovnicemi k efektivnímu řešení některých problémů. Kapitola 7 je věnována ilustrativním příkladům analytického a numerického řešení vybraných typů rovnic. Tato kniha rozšiřuje teoretický základ a poskytuje širší rozhled o souvislostech mezi různými oblastmi aplikace modelování a numerických simulací. Pro úplnost přidáni Maxwellové [9, 10, 8].

10

1 Úvod Na matematiku je možné pohlížet jako na soubor výrazových prostředků, druh jazyka, jehož prostřednictvím můžeme popisovat a předávat informace o světě, který nás obklopuje. Matematicky zaznamenané informace mají často formu rovnic. Tímto způsobem sestavený popis označujeme jako matematický model. Matematickým modelem se obvykle, z praktických důvodů snažíme o takovou prezentaci reality, která postihne klíčové vlastnosti modelovaného jevu, tj. jen ty, které jsou z hlediska zkoumání důležité. Dílčí děje s malým vlivem na vývoj a výsledný stav zkoumaného systému jsou při popisu zpravidla zanedbávány. Model ve výše uvedeném smyslu slova se skládá minimálně ze tří částí: 1. geometrie, 2. fyzikálních principů, 3. parametrů. Geometrie popisuje obecně prostorové rozložení všech součástí modelu, zejména jejich hranice. Fyzikální principy pak určují, jaké rovnice pro jednotlivé součásti modelu se použijí (Maxwellovy rovnice, zákony přenosu, principy mechaniky). Parametry popisují fyzikální vlastnosti jednotlivých částí modelu (permitivita, tepelná vodivost, modul pružnosti). Termínem modelování budeme v našem kontextu rozumět zjišťování fyzikálních vlastností reálných objektů a dějů jednak na základě vhodně zjednodušeného matematického popisu zkoumaných objektů a jednak aplikací zásadních fyzikálních principů. Stupeň zjednodušení popisu objektu i výběr fyzikálních principů závisí na tom, jaké informace chceme získat a s jakou jejich věrohodností se spokojíme. Typickou úlohou modelování je zjistit odezvu na podnět. Jednoduchým příkladem podnětu je (obecně prostorové) rozložení proudů, odezvou je pak buzené magnetické pole v prostoru. Jiným příkladem podnětu může být prostorového rozložení teploty. Odezvou je pak tepelný tok v prostoru. K popisu tepelného toku ovšem musíme znát i vlastnosti prostředí, které se obecně nazývají parametry úlohy. O tom jsme se zmínili již výše. V tomto případě je parametrem tepelná vodivost, která je obecně funkcí souřadnic. Z jednoduchých příkladů je zřejmé, že zadané veličiny, tj. podněty, i hledané veličiny, tj. odezvy, chápeme jako funkce souřadnic x, y, z, případně i času t, zapsáno matematicky u(x, y, z, t). Protože je veličina rozložena v prostoru, mluvíme o poli. Podle závislosti či nezávislosti veličiny na čase rozlišujeme stacionární nebo nestacionární procesy. Vývoj veličin u nestacionárních procesů závisí na čase. Zkoumané jevy jsou zpravidla popsány dvěma matematickými prostředky: 1. integrálními rovnicemi, 2. parciálními diferenciálními rovnicemi. Oba popisy v zásadě umožní zjistit odezvu na zadaný podnět. Výsledek je v obou případech stejný, pokud se jej podaří vyjádřit analyticky. Mezi oběma přístupy je ale zásadní rozdíl. Integrální přístup zahrnuje geometrii modelu v mezích integrálu. Diferenciální rovnice vyžaduje geometrii modelu v okrajových podmínkách. Bez ohledu na použitý postup platí, že obecné přesné řešení úlohy ve formě tzv. analytického řešení (tj. explicitní zápis hledané funkce u(x, y, z, t) matematickým vzorcem) je možné nalézt jen ve velmi jednoduchých speciálních případech, např. pravidelná a velmi jednoduchá tělesa, přičemž alespoň jeden rozměr je nekonečný (poloprostor, vrstva, nekonečný válec kruhového průřezu), nebo naopak alespoň jeden rozměr je nulový (bod, přímka, úsečka, kruh, rovina). Tyto modely jsou velmi omezující a obvykle nezahrnují dostatečně základní rysy zkoumaného objektu. V současné praxi nejběžnější a uznávaný postup je tzv. numerické řešení, tedy řešení přibližné, které spočívá v tom, že původní kontinuální (spojitá) úloha s nekonečně malými elementy 11

je nahrazena úlohou diskrétní s konečnými elementy. V případě integrálu jde o numerickou integraci, při diferenciálním přístupu se používá metoda konečných prvků (MKP či FEM - Finite Element Method). Ojediněle se můžeme setkat s tvrzením, že běžně používané rovnice jako jsou rovnice difuze, vedení tepla, elektromagnetismu apod. jsou jen přiblížením reality, mimo jiné už z důvodu reprezentace hmoty jako kontinua. O správnosti těchto rovnic svědčí především skutečnost, že byly ověřeny na základě experimentů a odvozeny z naměřených dat. Neexistuje jediný experiment, který by některou z nich vyvracel. Výsledky získané jejich řešením jsou v souladu s výsledky měření. Pokud jde o argument, že obvyklý popis výše uvedených dějů (difuze a další) nerespektuje diskrétní strukturu hmoty, pak je nutno si uvědomit, že důsledky diskrétní struktury se běžně nemohou projevit. Pro diskrétní strukturu hmoty platí principy kvantové fyziky. V nerelativistické kvantové mechanice je to například Schrödingerova rovnice. Pokud bychom chtěli zohledňovat diskrétní strukturu hmoty, museli bychom tuto rovnici řešit pro každou částici zvlášť. Poněvadž i v zrnku písku je nepředstavitelný počet částic, nestačily by k numerickému vyřešení této soustavy rovnic všechny spražené počítače světa a to ani ve vzdálené budoucnosti. Navíc by řešení nemělo valný význam, poněvadž by se vztahovalo k jednotlivým částicím. Měřitelný výsledek bychom získali průměrováním, což umí s nesrovnatelně nižšími náklady statistická fyzika. Nejjednodušším případem k aplikaci statistické fyziky je zředěný plyn. Rovnice popisující modelovaný děj jsou přiblížením k realitě v tom smyslu, že přesně neznáme parametry v nich vystupující. Může to být např. koeficient vedení tepla, či relativní permitivita. Obvykle je nahrazujeme konstantami a předpokládáme, že prostředí je homogenní a lineární. Ve skutečnosti jsou v prostředí odchylky jak od homogenity, tak od linearity, což se může projevit nesouhlasem mezi výpočtem a experimentem, pokud oba byly provedeny s vysokou přesností. To je však chyba modelu a nezpochybňuje použité zákony. Zpřesněním modelu lze dosáhnout lepší shody, ale s větším úsilím. Příprava modelu a vlastní modelování je náročná tvůrčí práce, která předpokládá detailní pochopení řešeného problému. V zásadě platí: čím složitější model použijeme, tím více informací získáme a tím bude i jejich přesnost vyšší. V praxi se však vynaložená námaha se složitějším modelem nemusí vrátit, pokud do detailu nerozumíme řešenému problému. Můžeme totiž přidat nepodstatný rys modelu. Úvodem do této problematiky se zabývá tato práce. V předchozích částech jsme několikrát použili termín fyzikální princip. Pokusíme se o jednoduché vysvětlení. Fyzikální princip je tvrzení, které se přímo nedokazuje, ale jehož pravdivost se potvrzuje tím, že všechny z něho logicky odvozené závěry jsou v souladu se skutečností. Fyzikální princip se popisuje v matematické formulaci, někdy se nepoužívají standardní matematické pojmy. Naproti tomu fyzikální zákon vychází z experimentálních dat, která popíše poměrně jednoduchou matematickou formou. V teoretické mechanice je zákonem zákon síly, tj. druhý Newtonův zákon. Principy pak jsou např. d’Alembertův princip, Lagrangeovy rovnice nebo Hamiltonův princip. Ten poslední používá málo známý matematický pojem funkcionál. Jinými příklady fyzikálních principů mohou být termodynamické věty, případně Maxwellovy rovnice. Často se (nejen v mechanice) pravdivost principu potvrzuje právě tím, že se z něho odvodí odpovídající zákon.

12

2 Fyzikální veličiny a zákony V této části se soustředíme na základy fyziky. Obecně popíšeme fyzikální veličiny a to z řady hledisek. Zejména se soustředíme na systematický popis složitějších veličin. Přístup je většinou fyzikální, i když matematická stránka se nezanedbává.

2.1 Souřadná soustava Řešení úlohy začíná volbou souřadné soustavy. Tu v principu volíme podle symetrie modelu, jsou zde však další nezanedbatelná hlediska. V zásadě je možné volit soustavu přímočarou a křivočarou. Křivočará soustava se volí v případě, kdy symetrie úlohy vede ke zjednodušení rovnic. V praxi používané křivočaré soustavy jsou válcová a kulová. Zejména válcová soustava se s výhodou používá, pokud model vznikl rotací kolem určité osy. V technické aplikaci je těchto případů hodně. Rovnice popisující model se však musí převést do nové křivočaré soustavy. Mají složitější tvar a řeší se specifickými metodami. Proto se tímto případem nebudeme zabývat. Přímočará, kartézská, pravoúhlá souřadná soustava se používá v převážné většině případů. Její nespornou vlastností je jednoduchost. To je také důvod, proč se přímočaré soustavy s úhly odlišnými od pravých zásadně nepoužívají, i když by v principu sloužily stejně dobře. Osy kartézské soustavy značíme symboly X, Y, Z, případně X1 , X2 , X3 nebo x1 , x2 , x3 . První způsob je názorný, druhý umožňuje rychlejší matematický zápis. Často se používají jednotkové → − → − → − − → vektory ve směru souřadných os (viz obr. 1(a)). Označují se např. i , j , k (vektor i je ve →, − → − → − → směru souřadné osy X atd.), nebo − e 1 e2 , e3 , přičemž vektor e1 je ve směru souřadné osy X1 atd. Kartézská souřadná soustava může být pravotočivá (viz obr. 1(a)) a levotočivá (viz obr. 1(b)). U pravotočivé soustavy je pořadí os takové, že při otáčení pravotočivého šroubu od kladné poloosy X ke kladné poloose Y se tento šroub pohybuje ve směru osy Z vzhůru. Osou budeme dále rozumět její kladnou poloosu, to již v dalším nezdůrazňujeme. Totéž platí pro kombinace Y, Z, X nebo Z, X, Y , které jsou vytvořeny cyklicky z počáteční kombinace. V posledním případě se pravotočivý šroub při otáčení od osy Z k ose X pohybuje ve směru osy Y . U levotočivé soustavy se pravotočivý1 šroub při otáčení od osy X k ose Y pohybuje proti směru osy Z.

(a)

(b)

Obrázek 1: Pravotočivá (a) a levotočivá (b) souřadná soustava. Pravotočivá a levotočivá soustava jsou na obr. 1. Levotočivá soustava vznikne z pravotočivé 1

Levotočivý šroub je obtížné si představit, poněvadž se v praxi používá vyjímečně.

13

prohozením dvou os, což má formální analogii se změnou smyslu otáčení trojfázového elektromotoru prohozením dvou fází. Při kreslení souřadné soustavy by se mělo vždy překontrolovat, zda je pravotočivá2 . Volbě souřadné soustavy je nutno věnovat péči. Souřadná soustava by se měla volit vždy s ohledem na symetrii úlohy. Při její optimální volbě je popis objektů jednoduchý. Tím se zjednoduší i formulace rovnic, jejich řešení je snadnější, takže se ušetří mnoho práce. Výběr samotného druhu souřadné soustavy (pravoúhlá, sférická, . . .) je nedostačující. Lze například zvolit několik (nekonečné množství) pravoúhlých souřadných soustav, které se budou vzájemně lišit polohou počátků a orientací (natočením) souřadných os. Mezi souřadnicemi bodů v různých soustavách existují transformační vzorce, které umožňují přepočet souřadnic mezi jednotlivými soustavami. Pokud se soustavy liší jen posunutím středů, nové a staré osy jsou rovnoběžné, jsou vzorce velmi jednoduché. Pokud jsou souřadné soustavy navzájem pootočeny, pak je ke transformacisouřadnic nutno použít směrových kosinů. Touto problematikou se nebudeme podrobně zabývat. Je však důležitá při práci s tenzory. Speciální transformací je přechod od pravotočivé k levotočivé souřadné soustavě.

2.2 Fyzikální veličiny Veškeré exaktně formulované fyzikální zákony jsou v různě složitém matematickém tvaru, od jednoduchého výrazu, např. Ohmův zákon, přes diferenciální formu, např. Maxwelovy rovnice, až k výrazům s velmi speciálními symboly a operacemi, kterým rozumí jen hrstka odborníků. Fyzikální zákon se od svého matematického zápisu liší tím, že použité symboly mají přesný fyzikální význam (nazýváme je fyzikálními veličinami) a zákon sám popisuje určitý fyzikální jev buď staticky, nebo dynamicky. Popis projevů modelovaného děje v čase a/nebo v prostoru se získá řešením příslušné rovnice. Fyzikální veličiny označují symboly, často doplněnými o různé indexy. Na rozdíl od matematiky má každá veličina svůj název, rozměr a jednotku3 . Mezinárodní soustava jednotek (SI) obsahuje 7 základních fyzikálních jednotek: metr - m (délka), kilogram - kg (hmotnost), sekunda - s (čas), ampér -A (elektrický proud), Kelvin - K (termodynamická teplota), mol - mol (látkové množství) a candela - cd (svítivost). Rozměry jednotek všech ostatních veličin se dají vyjádřit s použitím jednotek základních. Fyzikálních veličin je velké množství a dělí se do různých skupin podle určitých kriterií. Známé rozdělení je na extenzivní a intenzivní. Hodnoty extenzivních veličin se dají sčítat. Extenzivními veličinami jsou například hmotnost, objem, náboj, tepelná kapacita, energie. Sčítání hodnot extenzivních veličin dává po fyzikální stránce smysl. Typickým příkladem intenzivní veličiny je teplota. Důležitými intenzivními veličinami jsou veličiny extenzivní vztažené na jednotku délky a její mocniny (liniová, plošná a objemová hustota látky, náboje), na jednotku hmotnosti (měrná tepelná kapacita), jednotku času (výkon). Intenzivní veličiny jsou funkcemi polohy, pro extenzivní veličiny to platit nemůže. Další časté dělení je na integrální a diferenciální veličiny. Integrální jsou definovány v určité oblasti prostoru, diferenciální v daném místě. Ke každé integrální veličině existuje veličina diferenciální. Mezi oběma veličinami existuje jednoznačný vztah. Jako příklad uvedeme integrální − → veličinu elektrické napětí U a diferenciální veličinu intenzitu elektrického pole E . Je-li známo → − → → elekrické pole E (− r ), kde − r je polohový vektor, lze napětí mezi dvěma body A a B o souřadnicích 2

Osy se podaří prohodit při troše nepozornosti opravdu snadno. Při prezentaci pak mohou vzniknout nepříjemné dotazy ve smyslu ”Jaký má význam použití levotočivé souřadné soustavy?”. 3 Rozměr je např. délka a jednotka m, nebo rozměr je hmotnost a jednotka kg, rozměr je proud a jednotka A atd. Termín rozměr se používá málo, např. při fyzikálních úvahach, častěji se pracuje přímo s jednotkami, např. při rozměrové kontrole odvozeného vztahu.

14

→ → daných polohovými vektory − r A, − r B určit integrálem − →

∫r B U=

− →→ − E (− r )d→ r.

− → rA

Obrácený převod je možný pouze pomocí elektrického potenciálu, elektrické napětí je rozdíl potenciálu mezi dvěma body. Z pouhé znalosti napětí nelze určit elektrické pole. Informací je příliš málo. To platí i v dalších případech, např. mezi proudem a prostorovou proudovou hustotou. Je nutno podotknout, že mezi kategoriemi extenzivní/intenzivní veličiny a integrální/diferenciální veličiny není zásadní rozdíl. První členění je fyzikální, druhé matematické. Ve fyzikálním zákonu vystupuje jedna velična jako podnět. Ta je dána a v nejjednodušší formě je na pravé straně vztahu. Druhá veličina je odezva, tu z fyzikálního zákona určujeme a v jednoduchém případě bývá vlevo. Jako primitivní příklad je Ohmův zákon, I=

− → − → i = γ E.

U , R

(1)

V integrálním tvaru (vlevo) je podnětem napětí U a odezvou proud I, v difrencíalním tvaru → − − → (vpravo) je podnětem intenzita elekrického pole E a odezvou je prostorová proudová hustota i . → − → − Všechny čtyři uvedené veličiny U, I, E , i nazveme aplikační. V běžných (diferenciálních) rovnicích se podnět a odezva musí specifikovat složitěji a vystupují i další podmínky. Druhý typ veličin jsou veličiny materiálové (vlastnosti materiálu). V případě Ohmova zákona je to elektrický odpor R pro uvedený integrální tvar zákona a elektrická vodivost γ pro tvar diferenciální. Materiálové veličiny se mohou objevit v integrální i diferenciální formě fyzikálních zákonů, viz Ohmův zákon. Základní fyzikální zákony mají diferenciální tvar. V takovém případě jsou materiálové veličiny obecně funkcí souřadnic.

2.3 Matematický popis Na předchozí fyzikální dělení navazuje matematické dělení veličin podle počtu složek, které potřebujeme k jejich popisu. Uvedeme stručné, i když ne zcela přesné, základní definice. Oblast je část roviny (dvourozměrná oblast) nebo prostoru (trojrozměrná oblast) na které sledujeme, případně modelem popisujeme fyzikální jev. Skalár je jednorozměrná veličina určená svou velikostí, tj. jedním číslem (např. teplota, energie, náboj). Skalární pole je funkce, která každému bodu oblasti přiřazuje hodnotu skaláru (např. rozložení teplot na území ČR). Vektor je veličina určená nejen svou velikostí, ale též směrem či smyslem v prostoru, tj. ve zvoleném souřadném systému je popsána dvěma složkami v rovině a třemi v prostoru (např. síla, rychlost, magnetická indukce, teplotní gradient) Vektorové pole je funkce, která každému bodu oblasti přiřazuje hodnotu vektoru (např. pole rychlostí proudící kapaliny). Tenzor druhého řádu4 je veličina, která popisuje vztah mezi dvěma vektory. Je obecně vyjádřen devíti složkami (32 = 9), aplikačními tenzory jsou elastické napětí a deformace. Materiálovými tenzory jsou například tenzor elektrické nebo hydraulické vodivosti. Základ 3 je při uvádění počtu složek zvolen proto, že žijeme v trojrozměrném prostoru. Pro rovinu by byl základ 2. Matematicky lze uvažovat prostor libovolné dimenze. Fyzikální význam to však nemá. 4

Podrobněji se tenzory zabýváme v části 5.

15

Tenzorové pole je funkce, která každému bodu oblasti přiřazuje hodnotu tenzoru (např. rozložení elastického napětí v prostoru). Tenzor vyššího řádu popisuje vztah mezi tenzory nižšího řádu. Tyto tenzory jsou jen materiálové, příkladem tenzoru třetího řádu je piezoelektrický tenzor s 33 = 27 složkami, elastické moduly nebo koeficienty jsou tenzory čtvrtého řádu s 34 = 81 složkami. Zobecníme-li definici tenzoru, pak je skalár tenzorem nultého řádu, poněvadž má 30 = 1 složku, a vektor je tenzorem prvního řádu o 31 = 3 složkách. Základem tenzorového počtu je popis složek tenzorů při transformaci souřadné soustavy, čímž se rozumí její pootočení kolem počátku.

2.4 Materiálové prostředí Většina fyzikálních zákonů předpokládá materiálové prostředí. Reálné prostředí si obvykle poněkud zjednodušujeme, idealizujeme. Základní idealizací je homogenní prostředí, které má všude stejné vlastnosti. Všechny parametry, popisující toto prostředí, jsou stejné ve všech bodech. Vlastnosti nehomogenního prostředí se liší od místa k místu. Zavádí se též heterogenní prostředí, které se skládá z několika části (oblastí) s hranicemi, na kterých se alespoň některé parametry mění skokem. Se skokovou změnou parametrů se můžeme setkat třeba v případě geologických zlomů provázených vzájemným vertikálním posunem vrstev různých hornin. Jednotlivé části heterogenního prostředí mohou být homogenní nebo nehomohgenní. Stupeň homogenity závisí na tom, které fyzikální veličiny uvažujeme. Pokud jde o skaláry (hustota, teplota, tlak), obvykle jsou homogenní. V případě vektorů (magnetizace, polarizace) lze homogenitu předpokládat jen u geometricky velmi jednoduchých oblastí, např. koule může být homogenně zmagnetizována vnějším magnetickým polem. U tenzorů druhého řádu (elastická deformace či elastické napětí, viz kapitoly 5.2.1 a 5.2.2 ) je to ještě složitější. Zde se homogenita pouze předpokládá pro zjednodušení úlohy. Homogenní prostředí může být izotropní nebo anizotropní. U izotropního prostředí jsou jeho vlastnosti, parametry, stejné ve všech směrech. Vlastnosti anizotropního protředí jsou různé v různých směrech. Příkladem anizotropního protředí je ideální krystal, který má pravidelné uspořádání atomů5 . Vlastnosti ideálního (ale i reálného) krystalu se mění podle směru, ale systematicky, tedy krystal patřící do určité skupiny (krystalografické třídy) má vždy stejnou anizotropii. Směrovou odlišnost vlastností zajišťuje pravidelné uspořádání atomů. Nehomogenní prostředí má různé vlastnosti v různých místech. V důsledku toho má také různé vlastnosti v různých směrech. Vlastnosti prostředí ale závisí na směru nahodile. Totéž nehomogenní prostředí může vykazovat různou anizotropii, proto v případě nehomgenního prostředí nemá smysl o anizotropii vůbec mluvit. Prostředí vyplňuje určitou oblast. Jednoduché příklady oblastí jsou prostor, poloprostor, vrstva, nekonečný válec, nekonečný hranol, rovina. Tuto oblast v prvním přiblížení považujeme za homogenní, izotropní nebo anizotropní. Další základní vlastnost prostředí je to, jak reaguje na podnět. Pokud je odezva úměrná podnětu, nazýváme takové prostředí lineární. V případě složitější závislosti odezvy na podnětu mluvíme o nelineárním prostředí. Každé materiálové prostředí je nelineární, nelinearita se však projevuje až při silném buzení, které v řadě případů nelze realizovat, buď z důvodu, že nemáme tak výkonné zdroje vyvolávající podnět, nebo, že tyto silné zdroje prostředí degradují. Poslední základní rozdělení je dělení prostředí na měkké a tvrdé. U měkkého prostředí při návratu podnětu k nulové hodnotě klesne odezva též na nulovou hodnotu. U tvrdého prostředí si odezva zachováná určitou hodnotu i při odeznění podnětu. Zachování odezvy po odeznění podnětu označujeme jako hysterezi. Typickým tvrdým prostředím jsou feromagnetické látky. Mezi linearitou a hysterezí (měkké a tvrdé prostředí) není přímá souvislot. V praxi je však prostředí tvrdé rovněž prostředím nelineárním. Linerání tvrdé prostředí je idealizací. 5 Reálný krystal má určité odchylky od pravidelného uspořádání atomů: chybějící atomy (vakance), nadbytečné atomy (intersticiály), poruchy pravidelnosti (dislokace, nebo tepelné kmity atomů) atd.

16

Nejsložitějším protředím je prostředí s hysterezí. Je nelineární a tvrdé, navíc však neexistuje jednoduchá funkční závislost mezi podnětem a odezvou. Typickým příkladem jsou feromagnetické látky, u nichž je odezva na podnět popsána hysterezní smyčkou, tedy graficky. Danému podnětu odpovídají obvykle dvě odezvy. Která z nich nastane, závisí na předchozím působení, např. zda narůstalo, či klesalo. Říkáme, že látka má paměť. Nejjednodušší případ, kterým se budemem dále zabývat, je prostředí homogenní, izotropní, měkké a lineární. Je to velmi užitečná a jednoduchá idealizace a obyčejně bývá v praxi splněna s postačující přesností. Další idelaizací je prostředí homogenní, izotropní a tvrdé. Na podnět tedy vůbec nereaguje. Tato idealizace má své omezení, je to např. první přiblížení pro popis feromagnetické látky.

17

18

3 Vektorová algebra Použití vektorů ve fyzice využívá zákonitosti, vztahy, vlastnosti a definice zavedené v matematice. V této kapitole je poskytnut souhrn základních matematických pojmů a vzorců, které nalézají uplatnění při popisu vybraných fyzikálních veličin a jevů. Vazba níže představovaného aparátu na konkrétní aplikace není v rámci této kapitoly demonstrována.

3.1 Vektorový prostor V matematice se pracuje s N rozměrným Eukleidovským prostorem E (N ) , v němž je každý bod určen N-ticí reálných čísel6 , souřadnic. Z Eukleidovského prostoru se odvozuje N rozměrný vektorový prostor. Vektor je definován rozdíly souřadnic bodů v Eukleidovském prostoru. Tyto rozdíly nazýváme složky vektoru. → Prvky vektorového prostoru značíme fyzikálně šipkou nad symbolem7 , např. − v . Tento vektor → − má složky (v1 , v2 , ..., vN ) = v . Ve vektorovém prostoru jsou definovány níže uvedené operace a speciální prvky: 1. Součet dvou vektorů je opět prvkem vektorového prostoru, tedy je to vektor. Symbolicky → → − → − − C = A + B, prakticky se sečtou jejich složky Ci = Ai + Bi , i = 1, 2, ...N . − → → − 2. Součin reálného čísla a vektoru je vektor, tedy opět prvek vektorového prostoru, B = a A, neboli Bi = a Ai , i = 1, 2, ...N . − → − − − → → − → → 3. Existuje nulový prvek 0 s vlastností A + 0 = A pro každý vektor A. Nulový vektor má − → všechny složky nulové, 0 = (0, 0, ..., 0). − → 4. Na základě předchozí podmínky lze ukázat, že ke každému vektoru A existuje vektor − → − → − → − → − → − → opačný A n , takže A+ A n = 0 . Pokud A = (a1 , a2 , ..., aN ), pak A n = (−a1 , −a2 , ..., −aN ). Ve vektorovém prostoru jsou tedy definovány operace sčítání vektorů a násobení vektoru skalárem. Není definováno násobení vektorů tak, aby jejich součin byl vektor8 . Neexistuje tedy jednotkový prvek9 a nelze nalézt inverzní vektor k danému vektoru. Poznámka. Reálná čísla tvoří vektorový prostor podle předchozí definice. Zde je definováno násobení, jednotkový prvek a inverzní prvek. Jedná se ale o speciální případ vektorového prostoru. Požadavky výše uvedené nevylučují, aby vektory v určitých prostorech měly další speciální vlastnosti. Naproti tomu množina racionálních čísel (zlomky) vektorový prostor netvoří, poněvadž není definováno násobení reálným číslem. Racionální čísla jsou podmnožinou čísel reálných (Q ⊂ R). Pro operace sčítání vektorů a pro násobení vektoru reálným číslem platí tyto známé zákony: − − → → • komutativnost: → u +− v =− v +→ u, − → → → − → • asociativnost: → u + (− v +− w) = (− u +→ v)+− w, → → • asociativnost pro násobení skalárem: a(b− v ) = (ab)− v, → → → → → − → • distributivnost: (a + b)− v = a− v + b− v , a(− u +→ v ) = a− u + a− v. 6

Často používají jiné termíny a symboly. V tištěné literatuře se používá i jiné značení, často je to tučné písmo. Symbol šipka je však názorný a proto se přidržíme psané formy užívané ve fyzice. 8 Skalární součin, definovaný později, je reálné číslo. Výsledek vektorového součinnu popsaného v samostatné kapitole lze sice interpretovat jako vektor, ale jedná se pouze o formu zjednodušeného zápisu antisymetrického tenzoru. − →− → − → → − 9 Vektor o složkách 1 = (1, 1, ..., 1) není jednotkovým vektorem, poněvadž součin A 1 = A není definován. Je to řadový vektor jako každý jiný. 7

19

Požadavky na komutativnost a asociativnost jsou spojeny zavedením linerání kombinace, která je součtem násobků vybraných vektorů. Libovolný vektor lze vyjádřit pomocí jiných vek→ → → torů − v 1, − v 2 , ..., − v m násobených reálnými čísly a1 , a2 , ..., am jako ∑ − → → u = ai − v i. m

(2)

i=1

Symbolem m se rozumí počet použitých vektorů a nijak nesouvisí s dimenzí N vektorového prostoru. Vzniká otázka, zda lze takto vyjádřit každý vektor. K tomu účelu se využívají tzv. lineárně − − − nezávislé vektory. Vektory → v 1, → v 2 , ..., → v m v (2) jsou nezávislé, pokud tento vztah dává nulový vektor jen v případě, že všechna reálná čísla a1 , a2 , ..., am jsou nulová. Lze ukázat, že v N rozměrném vektorovém prostoru existuje právě N nezávislých vektorů, jejichž lineární kombinací lze vyjádřit libovolný vektor. Tyto vektory tvoří bázi vektorového prostoru. Nabízí se otázka, jak tuto bázi najít v nejjednoušší formě. Odpověď na tuto otázku bude naznačena v kapitole 3.2.1.

3.2 Vektory ve fyzice Ve fyzice pracujeme v trojrozměrném prostoru E (3) . Bázi trojrozměrného vektorového prostoru tvoří tři ortonormální (kolmé, jednotkové) vektory ve směru souřadných os X1 , X2 , X3 , − → e 1 = (1, 0, 0),

→ − e 2 = (0, 1, 0),

− → e 3 = (0, 0, 1).

(3)

→ Vektor − v o složkách (v1 , v2 , v3 ) lze vyjádřit lineární kombinací (2) jednotkových vektorů báze (3) jako →+v − → − → − → v = (v1 , v2 , v3 ) = v1 − e (4) 1 2 e2 + v 3 e3 . Geometricky vektor znázorňujeme orientovanou úsečkou, z níž je na první pohled zřejmá velikost i směr vektoru. Zásadní rozdíl od vektorů zavedených v matematice je v tom, že vektory ve fyzice mají fyzikální význam, tj. mají přiřazenu jednotku. Naproti tomu vektory v matematice jsou pouze geometrické. Zatímco matematické vektory nemají pevně daný počáteční bod, fyzikální vektory jej obvykle mají10 , jsou tedy umístěné v prostoru, mívají definované působiště. Jednotce délky vektoru přísluší jednotka fyzikální veličiny vektorem reprezentované. Pro velikost vektoru můžeme podle vztahu (5) psát11 v u 3 u∑ − → u = |u| = t u2i . (5) i=1

3.2.1 Skalární součin → → Skalární součin je operace, která dvěma násobeným vektorům − u, − v přiřazuje výsledek v podobě reálného čísla. Zápis a definice skalárního součinu má obvykle jeden z následujících tvarů. ∑ → − − − → u ·→ v = (→ u,− v) = ui vi = u1 v1 + . . . + uN vN N

(6)

i=1

→ → → − První tvar zápisu, − u ·− v , se používá ve fyzice. Druhý tvar zápisu, (− u,→ v ), je spíše matematický. 10 11

Případně je počáteční bod omezen na určitou oblast. Velikost vektoru píšeme stručně bez šipky nad symbolem.

20

Fyzikální význam skalárního součinu vychází z geometrické interpretace vektorů podle obr. 2 → − a pro vektory − u,→ v se píše ve tvaru → → W =− u ·− v = uv cos α,

(7)

→ → → → kde u = |− u |, v = |− v | jsou po řadě velikosti vektorů − u a− v , α je úhel, který vektory svírají, viz obr. 2. Na znaménku tohoto úhlu nezáleží, bere se vždy kladně.

Obrázek 2: Geometrický význam skalárního součinu. → → Podle obr. 2 je skalární součin součinem velikosti vektoru − u a velikosti průmětu vektoru − v − → do směru vektoru u . − → Učebnicovým příkladem skalárního součinu je práce vykonaná konstantní silou F na přímé − → dráze o délce s a se směrovým (tečným) vektorem t . V uvažovaném případě se jedná o součin → − − → − průmětu F cos α síly F do směru vektoru → s = s · t , tedy W = s · F cos α. Skalární součin se ve fyzice vyskytuje dosti často. S jeho pomocí se vyjadřují například veličiny s významem energie. Pro skalární součin platí následující zákony: − → − → • komutativnost: → u ·→ v =− v ·− u, → − → − • asociativnost pro násobení skalárem: a(− u ·− v ) = (a→ u)·→ v, → → → → → → → • distributivnost: (− u +− v)·− w =− u ·− w +− v ·− w. Skalární součin vektoru se sebou samým pak definuje velikost vektoru (Euklidovská norma vektoru) v uN √ u∑ − → − → − → u2i (8) u = |u| = u · u = t i=1

Vektory nazveme ortogonální, jestliže jejich skalární součin (6) je nulový. Lze ukázat, že ve vektorovém prostoru existuje báze tvořená ortogonálními vektory. Každý z těchto vektorů ještě vynásobíme reálným číslem tak, aby jeho velikost (8) byla jednotková. Tím získáme ortonormální souřadný systém.

PŘÍKLAD 1 → − Spočítejte skalární součin zadaných vektorů − u, → v √ √ √ √ → − → u = (1/ 2,√1/ 2), − v√= (−1/√ 2, 1/√ 2) → − − u ·→ v = 1/ 2 · (−1)/ 2 + 1/ 2 · 1/ 2 = 0 − → Nulová velikost skalárního součinu ukazuje na vzájemnou kolmost vektorů → u, − v.

21

3.2.2 Vektorový součin − → − → Vektorový součin je operace, která ze dvou násobených vektorů A, B (z euklidovského prostoru − → − → E 3 ) vytvoří třetí vektor C. Existuje několik forem zápisu vektorového součinu. Vektor C, který → − → − je vektorovým součinem vektorů A a B v tomto pořadí, je popsán následujícími vztahy. → − → − − − → − → → C = A × B = [ A, B],

C1 = A2 B3 − A3 B2 C2 = A3 B1 − A1 B3 C3 = A1 B2 − A2 B1

(9)

− → → − → − − → První forma zápisu vektorového součinu A × B se běžně používá ve fyzice, druhý tvar [ A, B] je spíše matematický. Protože je výsledkem vektorového součinu nový vektor, musíme u něj určit všechny tři fyzikální atributy: velikost, směr a smysl12 . − → → − 1. Velikost vektorového součinu | A × B| odpovídá obsahu rovnoběžníku se stranami − → → − A = A , B = B a je dána vztahem C = AB sin α

(10)

kde α je úhel sevřený vektory, viz obr. 3. 2. Vektorový součin je kolmý k rovině tvořené vektory v něm vystupujícími. Tedy je kolmý k oběma vektorům, viz obr. 3. 3. Smysl vektorového součinu je dán pohybem pravotočivého šroubu, otáčíme-li jím od prvního vektoru k druhému.

AxB

B

A B sin A

Obrázek 3: Fyzikální význam vektorového součinu. Vektorový součin je možné získat například jednou ze dvou podobně jednoduchých metod, které svým postupem odpovídají výpočtu determinantu matice o rozměru 3 × 3 prvky. Příslušná matice D je sestavená z trojice identifikátorů e1 , e2 , e3 směrů souřadných os na prvním a ze složek násobených vektorů na zbývajících dvou řádcích.   e1 e2 e3 (11) D =  A1 A2 A3  B1 B2 B3 12

Třetí vlastnost vektorového souučinu vyplývá z nekomutativnosti této operace. Budeme-li například provádět vektorový součin dvou vektorů ležících v rovině desky stolu, pak směr výsledku vektorového součinu bude kolmý k desce stolu. Smysl vektorového součinu vyplyne z pořadí násobených vektorů a určí orientaci výsledného vektoru nahoru nebo dolů.

22

Rozvoj determinantu matice podle řádku. Determinant matice D ze vztahu (11) lze vypočítat tzv. rozvojem determinantu podle řádku nebo sloupce. Matematický postup rozvoje podle prvního řádku je následující 3 ∑ det(D) = (−1)1+j · ej · det(D1,j ), j=1

kde symboly e1 , e2 , e3 identifikují směry jednotlivých bázových vektorů v používané souřadné soustavě. Výraz det(D1,j ) označuje determinant menší matice o rozměru 2 × 2, která vznikne z matice D vyškrtnutím 1. řádku a j-tého sloupce jak ukazuje Obrázek 4(a).

+ - + - + + - + (a)

(b)

Obrázek 4: Rozvoj determinantu podle prvního řádku.

Výpočet determinantu podle Sarrusova pravidla. Výpočtu determinantu v tomto případě předchází vytvoření tabulky, která vznikne opsáním prvních tří řádků matice D a jejich doplněním prvními dvěma řádky matice na místa 4. a 5. řádku tabulky. Po vytvoření tabulky je třeba vypočítat součiny těch jejích prvků, které leží v Obrázku 5 na vyznačených čarách. Součiny musí být následně vynásobeny mocninami čísla −1, které jsou indikovány znaménky zapsanými u začátků orientovaných úseček. Jednotlivé složky výsledného vektoru jsou tvořeny výrazy (resp. čísly), které představují koeficienty u identifikátorů směrů souřadných os e1 , e2 , e3 .

Obrázek 5: Schema pro výpočet determinantu.

Využití cyklického opakování indexů. Indexy v prvních třech symbolech ve vztazích (9) pro výpočet složek vektorového součinu se cyklicky opakují (1, 2, 3 → 2, 3, 1 → 3, 1, 2). Pokud zakryjeme na Obrázku 6 index té složky − → − → Ci , i ∈ {1, . . . , 3} vektorového součinu vektorového součinu A × B, kterou chceme počítat, pak zbývající dvě číslice procházené po směru hodinových ručiček identifikují složky násobených − → − → vektorů A, B v prvním členu na pravé straně hledaného předpisu. Druhý, záporný člen vzorce vznikne záměnou pořadí indexů ve členu prvním. Podobně jako u skalárního součinu posoudíme zákony jimiž se řídí vektorový součin: 23

1 3

2

Obrázek 6: Cyklické opakování indexů ve vektorovém součinu. → − − − → − → → • Nekomutativnost: Ze vztahů (9) plyne A × B = − B × A. Vektorový součin není komutativní. − → → − − − − → − → → → • Neasociativnost: Lze ukázat, že platí A × ( B × C) ̸= ( A × B) × C. 

   (B1 C3 − B3 C1 )A3 + (B1 C2 − B2 C1 )A2 −(A1 B3 − A3 B1 )C3 − (A1 B2 − A2 B1 )C2  (B2 C3 − B3 C2 )A3 − (B1 C2 − B2 C1 )A1  ̸=  −(A2 B3 − A3 B2 )C3 + (A1 B2 − A2 B1 )C1  −(B2 C3 − B3 C2 )A2 − (B1 C3 − B3 C1 )A1 (A2 B3 − A3 B2 )C2 + (A1 B3 − A3 B1 )C1

− → − − → − → − → → → → − − • Distributivnost: Platí ( A + B) × C = A × C + B × C stejně jako − → → → − − → − → − → − − → C × ( A + B) = C × A + C × B. Vektorový součin je obyčejně interpretován jako vektor. Je tomu tak proto, že má právě 3 unikátní, nezávislé složky. Ve skutečnosti jde ale o tenzor druhého řádu, tj. o matici s velice speciálními vlastnostmi. Tenzor je takzvaně antisymetrický. Označíme složky vektorového součinu symbolem se dvěma indexy Vij , i, j ∈ {1, . . . , 3}. V23 = A2 B3 − A3 B2 V31 = A3 B1 − A1 B3 V12 = A1 B2 − A2 B1

(12)

Pro výše neuvedené složky vektorového součinu platí Vij = 0 pro

i=j

Vij = −Vji

pro

i ̸= j

(13)

Prvky na diagonále uvažovaného tenzoru jsou rovny 0 a prvky symetrické podle diagonály mají stejnou velikost, ale opačná znaménka. V části 5.2.3 je antisymetričnost vektorového součinu demonstrována na příkladu malého mechanického otočení kolem bodu. Zásluhou antisymetričnosti má vektorový součin některé neobvyklé vlastnosti. Ve fyzice se tento druh vektorů označuje jako vektory axiální, aby se odlišil od běžných vektorů, které se nazývají radiální. Axiální vektory se vyskytují poměrně často, jsou jimi např. téměř všechny momenty, dále pak magnetizace, což je vlastně též moment magnetického dipólu.

3.3 Některé užitečné vztahy Základní vztahy mezi vektory již byly výše uvedeny. Důležité jsou zejména lineární kombinace a distributivní zákon. V praxi se ale používají některé složitější vztahy, které je dobré si pamatovat. Uvedeme dva nejdůležitější. Oba se týkají vektorového součinu. Pro kombinaci skalárního a vekorového součinu platí → − → − − → → → − → − − − − → → − → A · ( B × C) = C · ( A × B) = B · ( C × A)

(14)

→ − Vektory se cyklicky střídají. Geometrický význam je objem rovnoběžnostěnu o stranách | A|, − → → − | B| a | C|. Můžeme si jej představit jako zdeformovaný kvádr při zachování objemu, viz Obrázek 7. Podle části 5.2.1 jde o obecný prostý střih, což je kombinace čistého střihu a pootočení, obojí je ve třech rovinách. Příklad prostého střihu v jedné rovině je na obr. 8. Vektorové součiny 24

(a)

(b)

Obrázek 7: Fyzikální význam kombinace skalárního a vektorového součinu. v rovnostech (14) mají velikosti obsahů jednotlivých stěn a skalární součiny představují průměty vektorových součinů do směrů zbývajících vektorů. Druhý užitečný vztah se týká vektorového součinu tří vektorů − → − − → − → → → → → → → − − → − − → − − →− − A × B × C = A × ( B × C) = B( A · C) − C( A · B)

(15)

Vzorec (15) nemá jednoduché vysvětlení13 , ale používá se často. Platnost obou vztahů se dokazuje rozpisem, např. podle (9). Ověření je dosti náročné pokud jde o spotřebu papíru, ale jde k němu využít programy k provádění analytických výpočtů jako jsou placený Maple nebo volně dostupný Sage viz kapitolu 4.1. Potvrzení správnosti vztahu získané pomocí Sage ukazuje následující rovnice.  

   (B1 C3 − B3 C1 )A3 + (B1 C2 − B2 C1 )A2 , (A1 C1 + A2 C2 + A3 C3 )B1 − (A1 B1 + A2 B2 + A3 B3 )C1 , (B2 C3 − B3 C2 )A3 − (B1 C2 − B2 C1 )A1 ,  =  (A1 C1 + A2 C2 + A3 C3 )B2 − (A1 B1 + A2 B2 + A3 B3 )C2 ,  −(B2 C3 − B3 C2 )A2 − (B1 C3 − B3 C1 )A1 (A1 C1 + A2 C2 + A3 C3 )B3 − (A1 B1 + A2 B2 + A3 B3 )C3

Obrázek 8: Prostý střih (zkos) v rovině XY.

13 Mnemotechnicky se pamatuje dobře. Nalevo je abecední pořadí, na pravé straně je ”BAC”minus ”CAB”. BAC má lingvistický smysl a CAB se z něho získá přehozením krajních písmen.

25

26

4 Vektorová analýza Je-li hodnota skalární nebo vektorové veličiny funkcí souřadnic, pak mluvíme o skalárním nebo vektorovém poli. Ve fyzikální praxi to jsou funkce poměrně jednoduché, lze s nimi provádět jednoduché diferenciální a integrální operace. Jejich fyzikální aplikací se zabývají první části této kapitoly. Poslední podkapitola obsahuje čistě matematický přístup.

4.1 SageMath: snáze, rychleji, přesněji Pro lepší pochopení dále představovaných matematických nástrojů jsou následující kapitoly doplněny řešenými příklady. Analytické řešení i relativně jednoduchých matematických úloh může být velmi náročné na čas a pozornost. Energii vynaloženou na řešení rovnic popisujících vybraný jev lze ušetřit využitím některého z existujících softwareových nástrojů. V této a v dalších kapitolách je pro výpočty modelových příkladů využito volně šiřitelného programu Sage (resp. SageMath, viz [19]). SageMath je svobodný matematický software s otevřeným zdrojovým kódem a s licencí GPL. Je sestaven z mnoha existujících balíčků s otevřeným zdrojovým kódem: NumPy, SciPy, matplotlib, Sympy, Maxima, GAP, FLINT, R a mnoho dalších (viz [20]). Skrze programovací jazyk vycházející z Pythonu nabízí SageMath kombinaci silných stránek jednotlivých balíčků (viz [21]). Úkolem Sage je podle slov autorů: Vytvoření funkční, svobodné alternativy k programům ” Magma (viz [2]), Maple (viz [13]), Mathematica (viz [23]) a Matlab (viz [7]).“ Program Sage má bohatou online dokumentaci s řadou ukázek funkčních částí kódu. Ze stránek programu Sage je možné stáhnout jeho verze pro všechny tři nejrozšířenější operační systémy Windows, Mac OSX i Linux. Kromě toho jde k výpočtům v Sage využít také uživatelského rozhraní v prohlížečích internetových stránek.

4.2 Operátor nabla Ke zjednodušení a zpřehlednění v kapitole uváděných vzorců bude využíváno vektorového ope− → rátoru ∇, který je definován vztahem ( ) − → ∂ ∂ ∂ ∇= , , (16) ∂x ∂y ∂z Uvedený operátor je možné použít (aplikovat) různými způsoby na skalární nebo vektorové proměnné (veličiny). Aplikace operátoru nabla na různé typy proměnných přináší také různé → − druhy výsledků. Nejčastěji prováděné operace s využitím operátoru ∇ budou dále v kapitole předvedeny.

4.3 Skalární pole Typickým skalárem ve fyzice je potenciál φ, který je obecně funkcí souřadnic. Ve statickém případě lze tedy pro potenciál zapisovat jako φ(x, y, z), případně φ(x1 , x2 , x3 ). Skalární funkce závislé na prostorových souřadnicích definují takzvaná statická skalární pole14 , jinými slovy, definují rozložení skalární veličiny v prostoru. 4.3.1 Gradient Na skalární funkci, jejímž reprezentantem je zde potenciál φ(x, y, z), můžeme aplikovat operaci gradient, který je definován vztahem ( ) − → ∂φ ∂φ ∂φ G(x, y, z) = grad φ(x, y, z) = ∇φ(x, y, z) = , , . (17) ∂x ∂y ∂z 14

Potenciál může být i funkcí času, pak se jedná o dynamický případ. Tím se ale zde nezabýváme.

27

Pokud použijeme tenzorového značení souřadnic (kapitola 5), které je výhodnější pro kompaktní matematický zápis, vypadá definice gradientu následovně. ( ) − → ∂φ ∂φ ∂φ G(x1 , x2 , x3 ) = grad φ(x1 , x2 , x3 ) = ∇φ(x1 , x2 , x3 ) = , , (18) ∂x1 ∂x2 ∂x3 Operace gradient převádí skalární pole na vektorové. To je důvod, proč jsme v definicích (17) i (18) uvedli na klíčových místech souřadnice. V našem případě se ze skalárního potenciálového → − − → pole φ(x, y, z) vytváří pole vektorové G(x, y, z). Veličinu G nazveme intenzitou. Z definice (17) je zřejmý význam jednotlivých složek vektoru intenzity. Složky vektoru intenzity udávají rychlost změny potenciálu v daném bodě ve směru příslušné osy. To je též fyzikální význam složek intenzity. Např. v bodě o souřadnicích (xv , yv , zv ) je rychlost změny potenciálu ve směru osy Y rovna ∂φ ∂y = Gy (xv , yv , zv ). Z výše uvedeného je snadné odvodit rychlost Gu změny potenciálu ve směru libovolného jed→ − notkového vektoru − u = (ux , uy , uz ), u = ∥→ u ∥2 = 1 s počátkem v bodě (xv , yv , zv ). Rychlost Gu → − − změny potenciálu ve směru vektoru → u je průmětem intenzity G(xv , yv , zv ) do směru zvoleného vektoru. Za využití skalárního součinu (část 3.2), můžeme pro rychlost změny Gu potenciálu φ, → v daném směru vektoru − u psát vztah − → − → → → Gu = − u · G = u(· G cos(α) = − u (xv , yv) , zv ) · G(xv , yv , zv ) = → =− u · grad φ = u ∂φ + u ∂φ + u ∂φ , x ∂x

y ∂y

(19)

z ∂z

− → → kde α je úhel mezi jednotkovým vektorem − u a intenzitou G. Ze vztahu (19) je zřejmé, že pokud → je zvolený směr vektoru − u totožný se směrem gradientu potenciálu φ, pak je úhel α nulový, α = 0, takže cos(α) = 1. Z toho vyplývá, že rychlost změny potenciálu ve směru gradientu je největší. Největší je i intenzita. Odtud je zřejmý fyzikální význam gradientu potenciálu. Gradient skalární funkce udává směr ve kterém se skálární pole mění nejrychleji a velikost vektoru gradientu odpovídá rychlosti této změny. − Je-li jednotkový vektor → u kolmý k vektoru intenzity či gradientu, pak je úhel α = π/2, takže cos(α) = 0. Ze vztahu (19) je pak zřejmé, že rychlost změny potenciálu ve směru kolmém ke gradientu je v daném bodě nulová. Ze vztahů (17) i (18) je zřejmé, že gradient potenciálu a intenzita se nezmění, pokud k potenciálu přičteme konstantu. Potenciál tedy nemá fyzikální význam, poněvadž není určen jednoznačně. Je to však velmi užitečná matematická veličina. Obvykle se tato konstanta volí tak, že potenciálu přiřadíme nulovou hodnotu v nějakém významném místě v prostoru. V elektrostatice bývá tímto místem často nekonečno. Pak lze v trojrozměrném prostoru definovat plochy o konstantním potenciálu, které nazýváme ekvipotenciály. Ekvipotenciálu, které přísluší hodnota potenciálu φ0 , získáme řešením implicitní rovnice φ(x, y, z) = φ0 .

(20)

Ze vztahu (19) a z textu pod ním je zřejmé, že vektor kolmý ke gradientu udává směr, ve kterém se v daném bodu potenciál nemění. Tento vektor tedy musí ležet v tečné rovině k ekvipotenciále procházející daným bodem. Vektor k němu kolmý je podle vztahu (19) ve směru gradientu. Gradient potenciálu je tedy kolmý k ekvipotenciále v daném bodu. Z definice (17) lze v principu vypočítat potenciál, pokud je známa intenzita, neboli jeho gradient. To však ponecháne až do další části 4.4.1.

PŘÍKLAD 2 Příkladem potenciálových polí ve fyzice jsou pole gravitační a elektrostatické. Někdy se k nim připojuje i pole magnetostatické tvořené permanentními magnety. To lze však lépe popsat pomocí vektorového potencialu. Jako elementární příklad uvažujme elektrostatické pole tvořené 28

bodovým nábojem Q. Náboj je v počátku souřadné soustavy, takže pro jeho potenciál v bodu → o polohovém vektoru − r = (x, y, z) můžeme psát vztah Q Q 1 Q 1 1 → r)= = φ(− · = · √→ ·√ (21) − − → 2 4πϵ0 r 4πϵ0 4πϵ 0 x + y2 + z2 r · r √ → kde ϵ0 = 8, 85.10−12 F/m je permitivita vakua a r = |− r | = x2 + y 2 + z 2 je vzdálenost místa kde potenciál počítáme od bodového náboje Q. Potenciál je nepřímo úměrný vzdálenosti r od bodového náboje. Elektrostatický potenciál samotný nemá fyzikální význam. Fyzikální význam má však rozdíl potenciálu mezi dvěma body, což je elektrické napětí. Dále součin náboje a potenciálu je potenciální energie náboje v elektrostatickém poli, což je fyzikální veličina. Samotná potenciální energie závisí na volbě nulové hladiny, při fyzikálních úvahách jde však především o její změnu a ta je jednoznačná. Intenzita elektrického pole vyvolaného nábojem Q, který leží v počátku souřadné soustavy, je definována vztahem ( ) − − →→ Q x r y z Q → E (− r ) = − grad φ = . = 3 , 3 , 3 2 2 2 2 2 2 2 2 2 4πϵ0 (x + y + z ) 2 (x + y + z ) 2 (x + y + z ) 2 4πϵ0 r3 (22) Uvedené vyjádření intenzity se od obecné definice (17) liší jen znaménkem. Tato odlišnost je → zdůvodněna v části 4.4.5. Intenzita elektrického pole má směr polohového vektoru − r . Potenciál se tedy mění nejrychleji v tomto směru a rychlost jeho změny je nepřímo úměrná čtverci vzdálenosti − → r, E = | E | = k/r2 , kde k je konstanta. Ekvipotenciální plochy jsou soustředné koule. Intenzita má radiální směr, tj. směr poloměru kulových ploch.

4.4 Vektorové pole U vektorového pole je každá složka vektoru obecně funkcí tří souřadnic15 (x, y, z), případně − polohového vektoru → r = (x, y, z). To můžeme symbolicky zapsat např. takto → − − A(→ r ) = (Ax (x, y, z), Ay (x, y, z), Az (x, y, z)).

(23)

Složky vektorového pole bývají obvykle matematicky poměrně jednoduchými funkcemi souřadnic a proto je lze relativně snadno integrovat a derivovat. V trojrozměrném prostoru může být integrace vektoru prováděna podél křivky nebo plochy. V prvním případě mluvíme o účinku vektoru podél křivky, v druhém případě se jedná o tok vektoru plochou. Účinek vektoru podél křivky i tok vektoru plochou bude podrobněji vysvětlen v dalším textu. 4.4.1 Účinek vektoru podél křivky − →→ Nechť se v oblasti prostoru, kde je definováno vektorové pole A(− r ), nachází křivka (K). Body − → křivky (K) jsou určeny polohovým vektorem r = (x, y, z). Počáteční bod je dán polohovým → → vektorem − r 1 a koncový bod vektorem − r 2 . Orientace či smysl křivky jsou určeny pohybem od počátečního bodu ke koncovému. Ve fyzice má tato křivka jednoduché matematické vlastnosti, → − hlavně všude existuje jednotkový, tečný vektor t , tj. křivka má první derivaci podle všech → proměnných. Dílčí přírůstek polohového vektoru označme d− r = (dx, dy, dz). Pro jeho velikost můžeme psát vztah √ → dr = |d− r | = (dx)2 + (dy)2 + (dz)2 . (24) − → Jeho směr splývá se směrem tečného vektoru t křivky v daném místě. 15

Obdobně jako u skalárního pole vynecháváme případ, kdy je vektorové pole též funkcí času.

29

→ − → → V bodě křivky o polohovém vektoru − r lze pro velikost průmětu vektoru A(− r ) do směru − → tečného vektoru t použít mj. tyto vztahy − → → − At = A cos(α) = A · t ,

(25)

→ − − → − → − → kde A = | A| je velikost vektoru A, α je úhel mezi vektorem A a tečným vektorem t (viz obr. 9). − → Dílčí učinek vektoru A na elementu křivky (K) definujeme následujícími, ekvivalentními vztahy − → → − → − → At dr = A cos(α)dr = A · d− r = A · t dr (26) Tečkou je ve vzorci (26) označena operace skalárního součinu vektorů.

A t

At

r r2

r1

Obrázek 9: Dílčí účinek vektoru. → − − → → Účinek vektoru A(→ r ) podél křivky (K) s počátečním bodem − r 1 a koncovým bodm − r 2 pak definujeme vztahem − → ∫r 2 ∫ − → → − W = A·dr = At dr. (27) − → r1

(K)

Abychom zdůraznili, že jde o integraci podél obecné křivky (K), dáváme její symbol do závorek (K). Podobně později použijeme pro plochu symbol (S) a pro objem symbolem (V ). Nejde tedy o integrační meze. − →→ Pokud představuje vektor A(− r ) silu, pak je vzorec (27) známou definicí práce. Další, neobvyklé, členění kapitoly 4.4.1 je provedeno za účelem objasnění významů několika často užívaných označení fyzikálních polí. Nulový účinek vektoru podél křivky → − Nejprve se podívejme na možnosti kdy může integrál (28) vektoru A podél křivky (K) nabýt nulové hodnoty. ∫ At dr = 0. (28) W = (K)

30

Účinek vektoru na křivce může být nulový například pokud: → − 1. Vektor A má všechny složky nulové. 2. Délka křivky je nulová. 3. Vektor je v oblasti integrce všude kolmý na křivku. 4. Délka křikvy je nenulová. Křivka je uzavřená. Fyzikální pole je definováno rozložením ska→ − − lárního potenciálu φ a integrovaný vektor je gradientem daného potenciálu A = ∇φ(→ r ). Pole v takovém případě označujeme jako pole potenciálové. − →→ Neuvažujme první tři triviální případy a soustřeďme se na případ kdy vektor A(− r ) je − intenzitou, tedy gradientem potenciálu φ(→ r ), tj. − →− − A(→ r ) = grad φ(→ r)=

(

∂φ ∂φ ∂φ , , ∂x ∂y ∂z

) ,

(29)

− → viz definiční vztah (17). Pro přírůstek d→ r polohového vektoru − r můžeme psát − d→ r = (dx, dy, dz).

(30)

− → → Pro dílčí účinek A · d− r ve vztahu (27) dostaneme kombinací vztahů (29) a (30) výsledek ( ) − → − ∂φ ∂φ ∂φ → − → → − A · d r = grad φ( r ) · d r = dx + dy + dz = dφ. (31) ∂x ∂y ∂z → − Dílčí účinek vektoru A je tedy roven změně potenciálu dφ. Po dosazení do definičního vztahu (27) pak dostaneme − → − → ∫r 2 ∫r 2 → → W = At dr = dφ = φ(− r 2 ) − φ(− r 1 ). (32) − → r1

− → r1

Docházíme k důležitému závěru. Pokud je fyzikální pole popsáno potenciálem, pak účinek vektoru intenzity pole na dráze mezi dvěma body závisí jen na poloze těchto bodů a vůbec nezávisí na zvolené dráze, viz obr. 10.

Obrázek 10: Účinek vektoru na dráze mezi dvěma body v potenciálovém poli. Jestliže počáteční a koncový bod křivky v potenciálovém poli splývají, pak je účinek gradientu potenciálu na této křivce nulový. Pole s takovými vlastnostmi nazýváme pole konzervativní. 31

Například v konzervativním, zemském gravitačním poli, které je potenciálové, se na uzavřené dráze žádná energie (potenciální) nezíská ani neztratí. − →− Pomocí vztahu (32) můžeme ze známé intenzity pole G(→ r ) určit potenciál následujícím způsobem: → 1. Bod o polohovém vektoru − r 1 zvolíme takovým způsobem aby v něm byl potenciál φ → − nulový, tj. φ( r 1 ) = 0. Poněvadž jde o významný bod, označíme jeho polohový vektor jako − → r 0. → → 2. Změníme označení bodu s polohovým vektorem − r 2 (horní mez integrálu) na − r. 3. Po tomto přeznačení dostaneme vztah pro výpočet potenciálu v obecném bodě určeném → polohovým vektorem − r − → ∫r − → − → φ(− r)= Gd→ r. (33) − → r0

PŘÍKLAD 3 Jako velmi jednoduchý příklad uvažujme opět bodový náboj Q, který se nachází v počátku souřadné soustavy. V části 4.3 jsme pro jeho potenciál a intenzitu elektrického pole odvodili po řadě vztahy (21) a (22) → − →− Q 1 Q − r − φ(→ r)= , E (→ r)= . (34) 4πϵ0 r 4πϵ0 r3 Podle vztahu (32) dostaneme pro změnu potenciální energie malého bodového náboje16 q > 0 → → při přesunu z místa o polohovém vektoru − r 1 do místa o polohovém vektoru − r 2 vztah − → ( ) ∫r 2 − →→ Q 1 1 − − → − → → − ∆W = q E ( r )d r = qφ( r 2 ) − qφ( r 1 ) = q − . 4πϵ0 r2 r1

(35)

− → r1

→ → kde r1 = |− r 1 | a r2 = |− r 2 | jsou po řadě vzdálesnosti počátečního a koncového bodu od počátku. Pokud se zkušební náboj vzdaluje, r2 > r1 , je ∆W < 0. Energie soustavy klesá. Elektrická síla koná práci při vzdalování náboje q. Při přibližování náboje je r2 < r1 , ∆W > 0, energie soustavy vzrůstá. Práci na přibližování vykonává vnější síla.

PŘÍKLAD 4 Pokud známe intenzitu elektrického pole (34), můžeme podle vztahu (33) určit potenciál. Jen je třeba dát pozor na znaménko, viz vztah (22). Protože na integrační dráze v potenciálovém poli nezáleží, volíme tu nejjednoušší, polopřímku s počátkem v místě kde se nachází náboj Q, který pole budí. Dolní mez integrace je v nekonečnu, r0 → ∞. Pak můžeme vztah (33) postupně upravovat − →

− →

r ∫r − → − Q ∫ → φ(− r)=− E d→ r = − 4πϵ 0 Q = − 4πϵ 0 16

∫r ∞

− → r0

1 dr r2

− → r0

=

− → r − d→ r r3

Q = − 4πϵ 0

∫r r0

→ 1− r r3

r → − Q ∫ · t dr = − 4πϵ 0 r0

r dr r3

= (36)

Q 1 4πϵ0 r .

Tento náboj q nesmí ovlivnit pole buzené původním nábojem Q > 0, tj. musí být splněna podmínka q << Q.

32

− → → Vektory − r a t jsou v tomto případě rovnoběžné, proto můžeme v jejich skalárním součinu psát přímo jejich hodnoty. Úpravami vztahů pro výpočet potenciálu jsme dospěli ke známému výrazu, který je shodný s dříve uvedeným vzorcem (21).

Nenulový účinek vektoru podél křivky V obecném případě vektorového pole může být účinek vektorů podél uzavřené křivky různý od nuly − → ∫r2 At dr ̸= 0. (37) − r→ 1

Pro zajímavost poznamenejme, že pokud by tento vztah (37) neplatil, tj. práce po uzavřené dráze by byla vždy nulová, nemohly by existovat tepelné stroje, které konají práci cyklicky. Nebyly by tedy ani dnešní automobily. Na druhé straně je ale nutno dodat, že se zde jedná o další (termodynamické) zobecnění vztahu (37). Jedná se o stavový prostor a ne o prostor geometrický. Vztahy (28) a (37) mají názorné vysvětlení, pokud v nich jako vektorové pole uvažujeme − →→ → → → vektor rychlosti − v proudící kapaliny, tj. A(− r)=− v (− r ). Jestliže kapalina proudí jedním směrem (laminární proudění) je účinek rychlosti na libobovolné, uzavřené křivce nulový, tj. platí vztah (28). Ke každému kladnému příspěvku od souhlasných vektorů (vektor rychlosti a jednotkový tečný vektor ke křivce) lze nalézt v absolutní hodnotě stejně velký záporný přírůstek s opačnou orientací. Pokud však v kapalině vznikají víry, lze křivku volit např. podél dráhy některé částice ve víru a ve směru rotace. Pak je smysl obou vektorů souhlasný a platí vztah (37) s kladnou hodnotou integrálu. Takové vektorové pole nazveme pole vírové.

PŘÍKLAD 5 Dvourozměrné vektorové pole je popsáno vztahem − → A = (−y, x). Určete jeho cirkulaci (účinek podél křivky) na kružnici (K) s poloměrem R a se středem v počátku. Vysvětlení: Úkolem je vypočítat křivkový integrál ∫ w=

− →− → A ds =

(K)

∫

− → − → A · t ds.

(K)

Pro lepší pochopení zadání pomůže obrázek 11. → − • Kružnice (K) má v každém bodě tečný vektor t . Tečný vektor je podle konvence orientovaný do směru oběhu proti směru pohybu hodinových ručiček. − → • V každém bodě [x, y] na kružnici (K) má podle zadání počátek vektor A = (−y, x), který → − → je kolmý k polohovému vektoru − r = (x, y) a rovnoběžný s tečným vektorem t .

33

Obrázek 11: Cirkulace vektoru. Integraci je možné provést několika způsoby. − → → − 1. Jde využít rovnoběžnosti vektorů A a t spolu se znalostí obvodu kružnice. − → − → 2. Za pomoci transformace souřadnic vektorů A a t můžeme nahradit křivkový integrál obyčejným určitým integrálem podle jedné proměnné. Řešení 1.: − → → − Skalární součin vektorů A a t je díky jejich rovnoběžnosti a jednotkové velikosti tečného → − → − vektoru t roven velikosti A vektoru A, tedy − → − → A · t = A · t · cos(0) = A · t = A. − → Velikost vektoru A je na kružnici (K) rovna jejímu poloměru R. √ − → A = | A| = (−y)2 + x2 = R S touto znalostí jde křivkový integrál w zjednodušit následujícím způsobem. ∫ ∫ ∫ − → − → w= A · t ds = Rds = R ds (K)

(K)

(K)

Po vytknutí konstantní ∫ hodnoty poloměru před integrál zbývá vypočítat obvod kružnice s poloměrem R, protože ds = 2πR. (K)

− → Účinek vektorového pole A na kružnici s poloměrem R je 2 · πR2 . Řešení 2.: → − − → Souřadnice vektorů A a t s počátky na kružnici (K) je možné vyjádřit v polární soustavě → souřadnic za pomoci poloměru R a úhlu α sevřeného osou x a polohovým vektorem − r . K tomuto poslouží několik následujících vztahů. x = R · cos α y = R · sinα − → − → A = (−y, x) t = (− sin α, cos α) 34

→ − Transformaci souřadnic jednotkového vektoru t vyjasňuje následující obrázek.

Obrázek 12: Zavedení polárních souřadnic. Po transformaci souřadnic, před samotnou integrací podle α, je ještě potřeba správným způsobem nahradit ds výrazem s dα. K tomu jde použít vztah mezi délkou kruhového oblouku a úhlem vyjádřeným v obloukové míře. dα =

ds , ds = R · dα R

Pak už je samotná integrace snadná. ∫ w=

− → − → A · t ds = R

(K)

∫2π

− → → − A · t dα = R · R

0

∫2π (− sin α)2 + (cos α)2 dα = 2πR2 0

Získaný výsledek je pochopitelně stejný jako v případě využití prvního způsobu řešení. Pojmy nevírové a vírové pole používáme/zobecňujeme pro libovolný typ vektorů. Typickým zobecněním je intenzita elektrického pole. Statické elektrické pole je nevírové, zatímco dynamické (časově proměnné) elektrické pole obsahuje víry. Práce vykonaná vírovým elektrickým polem, které pohybuje nábojem po uzavřené křivce ve směru intenzity pole, je při každém oběhu kladná. Toho se využívá například ve speciálním a velmi jednoduchém urychlovačí částic, betatronu. 4.4.2 Tok vektoru plochou Přepokládáme poměrně jednoduchou plochu aplikovatelnou na fyzikální problémy. To mimo jiné znamená, že plocha má dvě strany. Vektorová analýza vyžaduje, aby plocha byla orientovaná, tj. každý její element musí mít nejen velikost, ale i směr a smysl. Proto každému bodu plochy − přiřadíme jednotkový normálový vektor → n . Pokud jde o jeho smysl, je několik možností. Je-li plocha uzavřená, míří normálový vektor ven z plochy. Obecně lze u plochy volit stranu 1 a stranu 2. V takovém případě míří jednotkový normálový vektor, podle konvence, ze strany 1 na stranu 2. Může být však definován smysl oběhu po okrajové křivce plochy. Pak směr normály souvisí se směrem oběhu po okrajové křivce pravidlem pravotočivého šroubu (pravé ruky). Uvedené zásady respektuje i definice orientovaného elementu plochy podle obr. 13. Pro úplnost zásady zopakujeme. Element plochy jako vektor je určen velikostí plochy a jednotkovým → normálovým vektorem − n. 1. Pokud není definován směr oběhu po okrajové křivce, zvolíme stranu 1 a stranu 2 elementu plochy, pak smysl normálového vektoru je ze strany 1 na stranu 2, viz obr. 13(a). 35

2. Pokud je určen směr oběhu po okrajové křivce elementu plochy, je smysl normálového vektoru je určen pravidlem pratočivého šroubu, viz obr. 13(b).

(a)

(b)

Obrázek 13: Orientovaný element plochy: 13(a) Obecně. 13(b) Pro zadaný směr oběhu. Pro element plochy reprezentovaný vektorem můžeme psát vztahy − → → dS = − n dS,

− → dS = |d S |.

(38)

− → → kde dS je velikost elementu plochy, − n normála k elementu plochy a d S vektor elementu plochy. Klíčovým vztahem je tok vektoru plochou či jejím elementem, který nazýváme elementární tok. Slovo tok“ zde chápeme v zobecněně. V případě, že vektorové pole má význam rychlosti ” proudící nestlačitelné kapaliny, jedná se skutečně o tok. Přesnějí jde o průtok vyjádřený jako objem kapaliny proteklý plochou za jednotku času. Obvykle jde ale o zobecněné případy toku, např. tok intenzity elektrického pole. − → → − Před definicí elementárního toku, tj. toku vektoru A elementem plochy d S , musíme vyjádřit → − velikost průmětu vektoru A do směru normály elementu plochy, kterou označíme An a nazveme normálovou složkou. Tu lze vyjádřit např. těmito ekvivalentními vztahy → − → An = A cos(α) = − n · A,

(39)

− → − → − → → kde A = | A| je velikost vektoru A a α je úhel mezi vektorem A a normálou − n k elementu plochy, viz například obr. 14. − → → − Elementární tok dΦ vektoru A elementem plochy d S je pak součin normálové složky An vek− → toru A a elementu plochy dS. S využitím vzorců (38) a (39) lze vyjádřit základní ekvivalentní vztahy (viz též obr. 14) − → − → → → − dΦ = An dS = A cos(α)dS = A · d S = A · − n dS,

(40)

− → − → − → → kde A = | A| je velikost vektoru A, α je úhel mezi vektorem A a normálou − n , An = A cos(α) − → − → − → je je průmět vektoru A do směru normály n , d S je vektor elementu plochy a dS je velikost tohoto elementu. Tečkou je značený skalární součin vektorů. Tok Φ vektoru A plochou (S) pak získáme integrací při použití vztahů (40). Vybereme-li první možnost vyjádření dílčího toku dΦ, dostaneme vztah ∫ ∫ → − − → n · AdS. (41) An dS = Φ= (S)

(S)

36

Obrázek 14: Tok elementem plochy. Symbolem (S) je naznačeno, že integrace se provádí po ploše (S). Pokud se volí strany plochy bez nějakých fyzikálních důvodů, není znaménko toku určeno jednoznačně. Na štěstí je tento případ ve fyzice poměrně vzácný. Obvykle je definován směr oběhu po okrajové křivce plochy a pak je normála všech elementů definována jednoznačně. Tok vektoru plochou má názorný význam v případě proudění nestlačitelné kapaliny. Pak je − → → → vektorem rychlost − v jejích částic ( A = − v ) a tok je objem kapaliny proteklý plochou. Pojem tok zobecňujeme i na další vektorová pole, jak uvidíme vzápětí. Vzorec (41) je platný i v případě toku vektoru uzavřenou plochou (S) ohraničujucí část prostoru (objem). Znaménko toku ven z uzavřené plochy (ven z objemu) je kladné.

PŘÍKLAD 6 Jako velmi jednoduchý příklad uvažujme opět bodový náboj Q, který se nachází v počátku souřadné soustavy. V části 4.3 jsme pro intenzitu jeho elektrického pole odvodili vztah (22) − → − − r Q → E (→ r)= , 4πϵ0 r3

(42)

→ kde − r je polohový vektor. Pro výpočet toku intenzity elektrického pole zvolme jako uzavřenou plochu kouli o poloměru r0 se středem v místě bodového náboje Q, tj. v počátku souřadné soustavy. Jednotkový, − → normálový vektor → n k této kouli má v každém jejím bodu směr polohového vektoru − r , takže − → − → − → platí r = | r | = n · r . Na kouli má navíc polohový vektor konstantní velikost r = r0 . S využitím těchto informací můžeme vztah (42) po dosazení do (41) postupně upravovat S = 4πr2 ∫ → → → → Q ∫ r0 − Q ∫ − n ·− n r ·− n = dS = 4πϵ0 dS = Φ= En dS = 4πϵ0 r03 r03 dS = 8πr dr (S) (S) (S) (43) r0 Q 1 ∫ Q Q 1 2 = 4πϵ0 r2 8πr dr = 4πϵ0 r2 4πr0 = ϵ0 , 0

protože

∫

0

0

dS = 4πr02 je plocha povrchu koule o poloměru r0 . Odvodili jsme speciální a velmi

(S)

jednoduchý případ obecného Gausssova zákona z elektrostatiky.

4.4.3 Divergence Divergence je vektorová operace, která převádí vektor na skalár. Z tohoto hlediska je inverzní k operaci gradient. 37

Existují dvě základní definice vektorové operace divergence: diferenciální, která je matematickým vyjádřením vhodným pro praktické výpočty, a integrální, která vyjadřuje fyzikální význam operace. Obě dvě níže uvedené definice divergence jsou vzájemně ekvivalentní (bude vysvětleno). − → Zápis operace divergence vektoru A může mít dvě zavedené podoby − → − → div A = ∇ · A. První forma zápisu je více využívana fyziky. Druhá forma zápisu se hojně objevuje v matematických textech. Diferenciální definice divergence − →− Pokud je vektorové pole popsáno vektorem A(→ r ), který má první parciální derivace v určité oblasti, pak je divergence tohoto vektoru definována vztahem − → ∂Ax ∂Ay ∂Az div A = + + . ∂x ∂y ∂z

(44)

Integrální definice divergence Integrální definice divergence vychází z toku vektoru plochou. Nechť malý objem (V ) je obklopen − →→ plochou (S). V tomto objemu je definováno vektorové pole A(− r ). Definujme normovaný tok → − − → vektoru A( r ) uzavřenou plochou (S) vztahem ∫ An dS Φn =

(S)

(V )

,

(45)

− → → kde An je průmět vektoru A do směru vnější normály − n plochy (S) obklopující objem (V ). − → Tok Φ vektoru A uzavřenou plochou (S) je nenulový v případě, že převažuje výtok nad vtokem nebo naopak. Pokud převažuje výtok, říkáme, že v objemu se vyskytují zdřídla. Když převažuje vtok, hovoříme o výskytu nor17 . Pokud by vektorové pole bylo homogenní, jednalo by se o výtok vektoru z jednotkového objemu (nebo o vtok do něho). V případě nehomogenního vektorového pole můžeme mluvit o hustotě zdrojů, kterými jsou buď zřídla nebo nory18 . Integrální definici divergence dostaneme jako limitní případ definice (45), pokud objem (V ) a všechny jeho rozměry konvergují k nule ∫ An dS → − (S) , (46) div A = lim (V ) (V )→0 Z této definice je zřejmé, že dělení objemem je v definici divergence nutné. Jinak by limita (46) byla nulová, s výjimkou singularit19 . Integrální definice sice nevyžaduje existenci parciálních derivací vekrorového pole, musí však existovat limita (46). V běžné praxi jsou oba požadavky ekvivalentní, derivace je vlastně také limita. Z integrální definice (46) plyne fyzikální význam vektorové operace divergence. Popisuje normovaným způsobem výtok vektoru v daném bodě nebo vtok vektoru do něho. Divergence 17

Používají se také termíny propad nebo stoka. Všechny tyto pojmy jsou převzaty z mechaniky tekutin. Zde je zobecňujeme. 18 Zdrojem zde rozumíme jak zřídlo, tak noru, prostě případ, kdy je divergence nenulová. Zdroje rozlišujeme na zřídla nebo nory podle znaménka příslušejícího toku. Tato termionologie není ideální, poněvadž mezi zdrojem a zřídlem se obvykle nerozlišuje, Běžně znamenají jedno a totéž. 19 Singularitou je např. bod, či množina bodů, v nichž vektor v absolutní hodnotě roste nade všechny meze.

38

tedy udává hustotu zdrojů vektoru. Pokud je v nějaké oblasti divergence nulová, nejsou v této oblasti ani zřídla ani nory. Pak mluvíme o nezřídlovém poli. Pro lepší pochopení uvedeme několik velmi jednoduchých příkladů. V každém z příkladů uva→ − → žujeme proudění nestlačitelné kapaliny, kterému přísluší rychlostní pole − v (− r ), kde → r = (x, y, z) je polohový vektor. Pro jednoduchost v příkladech pracujeme s kulovou oblastí o objemu (V ) se středem v počátku souřadné soustavy. Nezabýváme se tím, zda jsou tyto přepoklady alespoň přibližně realizovatelné. Pokud je rychlost proudění konstantní, je divergence podle vztahu (44) nulová.

PŘÍKLAD 7 Pokud proud vychází z počátku souřadné soustavy, jeho rychlost vykazuje kulovou symetrii a narůstá lineárně, tj. − → − → v (→ r)=a− r = (a · x, a · y, a · z), (47) kde a je konstanta, je divergence konstantní → div(− v ) = 3a.

(48)

To znamená, že konstantní hustota zdrojů v oblasti vyvolá lineární závislost rychlosti proudění na vzdálenosti od počátku. Se zobecněnou aplikací tohoto případu se setkáváme např. v elektrostatice, kde platí Maxwellova rovnice − → div D = ρ, (49) − → kde D je elektrická indukce a ρ je objemová hustota volného náboje. Pokud je tato hustota náboje ve zvolené oblasti konstantní, narůstá v ní elektrická indukce lineárně. To je dost zvláštní případ, obvykle elektrické veličiny klesají20 .

PŘÍKLAD 8 Pokud bude √rychlost z předchozího zadání ubývat nepřímo úměrně vzdálenosti od počátku − r = |→ r | = x2 + y 2 + z 2 podle vztahu ( ) − → → − → v (− r ) = b rr2 = b x2 +yx2 +z 2 , x2 +yy2 +z 2 , x2 +yz2 +z 2 , (50) − → → v = |− v | = b| r | = b1, r2

r

→ kde b je konstanta, pak je divergence rychlosti − v dána vztahem b → div(− v) = 2. r

(51)

Z posledního vztahu plyne toto tvrzení: Pokud ve zvolené kulové oblasti mohutnost zdrojů klesá se čtvercem vzdálenosti od počátku, rychlost proudění nestlačitelné kapaliny klesá nepřímo úměrně této vzdálenosti.

Gaussova věta Pro vektorové pole s definovanými parciálními derivacemi v objemu (V ) obklopeném plochou (S) má významné postavení Gaussova věta21 , která vyjadřuje ekvivalenci diferenciální a integrální definice divergence (viz (44) a (46)). 20 21

To ovšem platí v tomto případě i ve velké vzdálenosti od zvolené oblasti Někdy se v literatuře uvádí též jako Gaussova-Ostrogradského věta.

39

Gaussova věta: − → − → Objemový integrál z divergence vektoru A je roven toku vektoru A uzavřenou plochou (S) ∫ ∫ − →− − − → → . ∇ · AdV = A→ n dS, (52) (V )

(S)

kde uzavřená plocha (S) je povrchem tělesa o objemu (V ). Formálně lze Gaussovu větu 52 přepsat do následující podoby ∫ ∫ − → div AdV = An dS. (V )

(53)

(S)

V případě proudění nestlačitelné kapaliny vyjadřuje Gaussova věta zákon zachování hmotnosti. Kapalina, která je v objemu produkována zřídly musí odtékat a kapalina pohlcovaná norami musí odněkud přitékat. Pokud mají integrály v Gaussově větě zápornou hodnotu, pak převládá vliv nor a převažuje vtok vektoru do oblasti. Naopak, když vychází integrály v Gaussově větě kladné, pak převládá vliv zřídel a převažuje výtok vektoru z oblasti. Pokud v objemu nejsou ani zřídla ani nory, je celkový tok okrajovou plochou nulový. Co do objemu vteče, musí také vytéci. Názorná ilustrace Gaussovy věty pro dvourozměrný případ je zachycena na obr. 15. Uvnitř objemu se uvažuje jedno zřídlo a dvě nory. Tok okrajovou plochou je určen blízkým zdrojem či norou.

Obrázek 15: Ilustrace tvrzení Gaussovy-Ostrogradského věty.

4.4.4 Rotace Rotace je vektorová operace, která převádí vektor na (ještě složitější) vektor. Ze tří vektorových operací, gradient, divergence a rotace, je nejsložitější. Stejně jako u divergence existují dvě základní definice vektorové operace rotace. Jsou to definice diferenciální a integrální. Diferenciální definice rotace → Uvažujme vektorové pole A(− r ), které má ve zkoumané (modelem popisované) oblasti definovány

40

→ první parciální derivace. Rotace vektoru pole A(− r ) je pak definována vztahem ( ) − → − → ∂Ay ∂Ax ∂Az ∂Ay ∂Az ∂Ax rot A = ∇ × A = − , − , − . ∂y ∂z ∂z ∂x ∂x ∂y

(54)

− → − → Jedná se o diferenciální definici. Rotace je vektorový součin operátoru ∇ a vektoru A. Pro jeho provedení lze využít schematu uvedeného v následujícím příkladu.

PŘÍKLAD 9 − → − → Určeme x-ovou složku Bx vektoru rotace B = ∇ × A. 1. Na níže uvedeném obrázku zakryjeme označení souřadnice vektorového součinu, kterou počítáme, tedy x.

x z

y

Obrázek 16: Schéma výpočtu rotace vektoru. 2. Postupujeme ve směru šipek, tj. proti směru hodinových ručiček, po nezakrytých značkách. 3. Pořadí procházených značek nám říká, že v uvažovaném případě budeme derivovat složku − → Az vektoru A podle proměnné y. 4. Proměnné y, z následně prohodíme. Provedeme parciální derivaci Ay podle z. − → 5. Výsledek, x-složku vektoru rotace ∇ × A získáme odečtením druhé, provedené parciální ∂Ay z derivace od první, tj. Bx = ∂A ∂y − ∂z . − → Analogicky k uvedenému postupu lze určit i zbylé dvě složky By , Bz vektoru B. Pozornému čtenáři neujde podezřelý směr procházení vrcholů trojúhelníku ve schematu. Ten − → je opačný než v případě počítání vektorového součinu. Vysvětlení je jednoduché. Vektor A stojí ve vektorovém součinu s operátorem ∇ na druhém místě. Vektorový součin není komutativní a obrácení směru procházení vrcholů zajistí získávání indexů v zavedém pořadí: 1) derivovaná složka vektoru, 2) nezávisle proměnná. Stejným způsobem jako v uvedeném příkladu lze vypočítat i ostatní složky rotace vektoru → − A s tím rozdílem, že ve využitém schematu zakryjeme příslušné písmeno. Integrální definice rotace Integrální definice rotace vychází z účinku vektoru na křivce. Nechť má malá oboustranná plocha (S) na okraji křivku (C), na které je stanoven směr oběhu. Na ploše (S) budiž definováno − →→ − →− vektorové pole A(− r ). Normovaná normálová složka rotace (neboli víru) vektoru A(→ r ) na uzavřené křivce (C) je předepsána vztahem ∫ At dr Vn =

(C)

(S)

41

,

(55)

→ − − → kde At je průmět vektoru A do směru tečny t křivky (C) na okraji plochy (S). Pokud je vír Vn 22 (C) orientován ve směru oběhu křivky, potom je integrál (55) kladný. V případě orientace víru proti směru křivky je integrál záporný. Pokud by vektorové pole bylo homogenní, jednalo by se ve vztahu (55) o popis víru pro jednotkovou plochu. V případě nehomogenního vektorového pole můžeme mluvit o hustotě vírů. Je vhodné si uvědomit, že integrální definice se netýká vektoru rotace, ale jen jeho složky Vn ve směru normály malé plochy ohraničené křivkou (C). Tato složka je skalární stejně jako účinek vektoru po hraniční křivce malé plochy. Integrální definice rotace je limitním případem vztahu (55). Plochu (S) a všechny její rozměry necháme konvergovat k nule. Tím získáme vzorec pro normálovou složku rotace ∫ At dr − → − → − (C) n · rot A = lim , (56) rotn A = → (S) (S)→0 → Vektor − n je jednotkový vektor normály k nekonečně malé ploše dS. Jeho smysl je dán směrem oběhu po okrajové křivce (C) elementu plochy dS podle pravidla pravotočivého šroubu. Zdůraz→ něme, že definice (56) vyjádřuje pouze průmět rotace do směru normály − n k elementární ploše dS. Ze vzorce (56) plyne fyzikální význam vektorové operace rotace. Popisuje normovaným způ− → sobem vír vektoru v daném bodě. Pokud je v nějaké oblasti rotace rotn A rovna nule, nejsou v této oblasti víry. V takovém případě mluvíme o nevírovém poli.

PŘÍKLAD 10 Jako jednoduchý příklad rotace uvažujme tenkou desku, která rotuje kolem osy z úhlovou rych− lostí ω. Vektor rotace má pak složky → ω = (0, 0, ω). Typickou technickou realizací je pevný disk v počítači nebo disk DVD. Pro obvodovou rychlost ve vzdálenosti r od počátku platí zvtah → → − → → → → v =− ω ×− r , v = ω · r, r = |− r |, v = |− v |, ω = |− ω |.

(57)

Výzam použitých značek osvětluje následující obrázek.

Obrázek 17: Ilustrace popisu rotačního pohybu DVD. → Vektor úhlové rychlosti − ω leží v počátku souřadné soustavy a vystupuje z papíru směrem vzhůru. 22

Tento termín byl opět jsou převzat z mechaniky tekutin a zde je zobecněn.

42

K výpočtu rotace potřebujeme složky obvodové rychlosti v obecném bodě, který je popsán buď kartézskými souřadnicemi (x, y), nebo polárními souřadnicemi (r, φ). Pro složky obvodové rychlosti odvodíme vztahy vx = −v sin φ, vy = v cos φ, vz = 0.

(58)

Goniometrické funkce ve vztahu (58) lze vyjádřit jako sin φ =

y , r

cos φ =

x . r

(59)

Po dosazení vzorců (57) a (59) do (58) dostaneme výsledné vztahy vx = −yω, vy = xω, vz = 0.

(60)

Odtud již snadno spočítáme složky rotace s využitím definičníchh vztahů (54) → rotx − v =

∂vy ∂vz − = 0, ∂y ∂z ∂vx ∂vz → roty − v = − = 0, ∂z ∂x ∂vy ∂vx → v = rotz − − = ω + ω = 2ω. ∂x ∂y

(61)

Rotace má pouze složku ve směru osy z, je konstantní na celém disku a co do velikosti je dvojnásobkem úhlové rychlosti. Rotace tedy velmi těsně souvisí s vektorem úhlové rychlosti. Příklad s rotujícím diskem je sice jednoduchý, není však typický pro vektorové pole vykazující rotaci. Konstantní rotace je zajištěna tuhostí disku. Typickým příkladem rotace je vír na vodní hladině. Pokud by rotace měla být konstantní, musel by se vír s rostoucí vzdáleností od středu otáčet stále vyšší obvodovou ryhlostí, což nepozorujeme. Rotace může být konstaní pouze v oblasti blízké jeho středu, pak klesá. Realističtější je v prvním přiblížení uvažovat, že obvodová rychlost klesá nepřímo úměrně vzdálenosti od středu víru. Při splnění tohoto předpokladu bude rotace klesat se čtvercem vzdálenosti, podobně jako tomu bylo u divergence. Konstantní rotace stejně jako konstantní divergence v případě rychlostního pole vyžaduje, aby rychlost narůstala lineárně se vzdáleností od počátku, viz vztah (60) pro rotaci a vztah (47) pro divergenci23 . Jediný rozdíl je v tom, že v případě divergence má rychlost směr poloměru, tedy radiální (směr normálového vektoru), zatímco u rotace má směr kolmý na poloměr, tj. axiální (směr tečného vektoru). To pro dvourozměrný případ názorně ukazuje obr. 18. Na obrázku 18(a) mají vektory pole radiální směr. Pole je zřídlové a nevírové, protože rychlost nemá axiální složku. Vektory s počátky na některé ze soustředných kružnic se středem v počátku mají stejnou velikost, která narůstá lineárně s poloměrem. Z toho důvodu je divergence vektorového pole konstantní (ve vyšetřované oblasti). 23

To by nemělo překvapit, poněvadž se jedná o tytéž derivace, jen v jiných vztazích.

43

(a)

(b)

(c)

Obrázek 18: Porovnání vektorového pole pro konstantní divergenci a rotaci. 18(a) Zřídlové a nevírové pole, 18(b) Vírové a nezřídlové pole, 18(c) Zřídlové i vírové pole. Na obrázku 18(b) je rychlost kolmá k poloměru a tedy tečná ke kružnicím se středem v počátku. Protože zobrazené vektory nemají radiální složky, pole je nezřídlové, ale je vírové. Rychlost má i v tomto případě konstantní velikost na soustředných kružnicích a roste úměrně jejich poloměru. Rotace je tedy konstantní ve vyšetřované oblasti. Obrázek 18(c) je kombinací obou předchozích. Vektory mají jak radiální, tak axiální složku. Vektorové pole je zřídlové i vírové. Rotace i divergence jsou konstantní, pprotože obě složky mají konstantní velikost na soustředných kružnicích a rostou spolu s poloměrem. Uvedené charakteristiky pole (zřídlovost, vírovost) pole lze posuzovat i v trojrozměrném prostoru pokud kružnice nahradíme soustřednými koulemi.

PŘÍKLAD 11 Zřídlové a nevírové elektrické pole vytvářejí elektrostatické náboje. Je to pole statické. Vírové elektrické pole vzniká pouze při časové změně magnetického pole. Je to pole dynamické. Tato asymetrie je způsobena skutečností, že neexistují magnetické náboje.

Stokesova věta Velmi důležitá je pro modelování fyzikálních polí Stokesova věta, která vyjadřuje ekvivalenci diferenciální a integrální definice rotace (viz (54) a (56)). Stokesova věta: − → − → Tok vektoru rot A plochou (S) je roven cirkulaci vektoru Apo uzavřené křivce (C), ∫ ∫ − → − →→ − → → − − . = ( ∇ × A) · n dS A t dl, (62) (C)

(S)

kde plocha (S) je ohraničena křivkou (C). Stokesovu větu (62) je možné formálně vyjádřit také jako ∫ ∫ − → At dr. rotn AdS = (C)

(S)

44

(63)

Plocha (S) má na okraji křivku (C). Smysl oběhu po této křivce je dán jednotkovým tečným − → → vektorem t . Symbolem − n je označen jednotkový normálový vektor plochy (S), jehož smysl je určen podle pravidla pravotočivého šroubu smšrem oběhu po křivce (C). Na levé straně integrálu (62) stojí součet normálových složek elementárních vírů na ploše (S). Na pravé straně je vír pro okrajovou křivku (C). K vysvětlení Sotesovy věty slouží obr. 19. Na ploše (S) jsou dva elementární pravotočivé víry a jeden vír levotočivý. Na obrázku je též naznačen vliv elementárních vírů na okrajový vír.

Obrázek 19: Ilustrace tvrzení Stokesovy věty. Vazba mezi tokem rotace vektoru plochou a účinkem vektoru na hranici plochy je demonstrována na následujícím příkladu.

PŘÍKLAD 12 Trojrozměrozměrné vektorové pole je popsáno vztahem → − A = (−y, x, 0). Určete cirkulaci (účinek podél křivky) vektorového pole na kružnici (K) ležící v rovině xy, mající poloměr R a střed v počátku souřadné soustavy. Řešení: − → Nejprve vypočítáme rotaci vektoru A. − → → − → − ( ) ( ) ( ) i j k → ∂ → − ∂x ∂y → − ∂y − → − − ∂x → ∂ ∂ ∇ × A = ∂x ∂y ∂z = i − + j + k + = (0, 0, 2) ∂z ∂z ∂x ∂y −y x 0 − → − → Výsledný vektor ∇ × A zintegrujeme na ploše (S) ohraničené kružnicí (K). Tato plocha má → jednotkový normálový vektor − n = (0, 0, 1). ∫ (S)

→ → − → − ( ∇ × A) · − n dS =

∫

S = πr2 (0, 0, 2) · (0, 0, 1)dS = dS = 2πrdr

∫R = 2 · 2π rdr = 2πR2 0

(S)

→ − → − Tok vektorového pole ∇ × A plochou ohraničenou kružnicí (K) vyšel rovný účinku vektoru ∫ → → − − → − A na hranici této plochy (viz dříve získaný výsledek A · t ds = 2πR2 ). (K)

45

4.4.5 Potenciály V předchozích částech jsme popsali tři základní operace vektorové analýzy: gradient (viz kapitolu 4.3.1), divergence (viz kapitolu 4.4.3) a rotace (viz kapitolu 4.4.4). Nyní shrneme jejich základní vlastnosti a naznačíme možnosti využití. Navíc pro případ skalárního pole zde odvodíme potenciál obecněji než v úvodní části 4.3 této kapitoly. Základní vlastnosti operací gradient, divergence a rotace jsou tyto 1. Gradient převádí skalární pole na vektorové. 2. Divergence transformuje vektorové pole na pole skalární. 3. Rotace vytváří z vektorového pole opět vektorové pole s jinými vlasnostmi. Z toho vyplývá, že ani operaci gradient, ani divergenci nelze opakovat. Např. kombinace grad(grad φ) aplikuje v druhém kroku operaci gradient na vektor a nikoliv skalár, což je chyba. Jedinou operací, kterou lze opakovat vícenásobně, je rotace. Aplikuje se na vektor a jejím výsledkem je opět vektor. Různý typ výstupu operací vektorové analýzy je též důvodem, proč jednotlivé operace nelze libovolně kombinovat. Přípustné jsou pouze ty kombinace, které uvádíme v systematickém pořadí: − → • grad(div A), • div(grad φ), − → • div(rot A) a • rot(grad φ). → − První z uvedených operací (grad(div A)) dává poměrně složitý vztah, který nemá významné fyzikální využití. Naproti tomu ostatní operace jsou velmi jednoduché, mají názorný význam (bude vysvětlen v kapitole 4.6) a nachází časté uplatnění ve fyzice. 1. Divergence gradientu vede na vztah ∇ · (∇φ) = div(grad φ) = kde symbolem

( ∆=

∂2φ ∂2φ ∂2φ + + 2 = ∆φ, ∂x2 ∂y 2 ∂z

∂2 ∂2 ∂2 + + ∂x2 ∂y 2 ∂z 2

(64)

) (65)

je označen známý Laplaceův operátor. 2. Divergence rotace je nulová, jak se snadno přesvědčíme dosazením → − → − ∇ · (∇ × A) = div(rot A) = 0.

(66)

3. Rotace gradientu je rovněž nulová. ∇ × (∇φ) = rot(grad φ) = 0.

46

(67)

Skalární potenciál Nulovost rotace gradientu skalární funkce φ indikuje nevírovost potenciálového pole. Následující − →− dvě rovnice (68) a (69) umožňují zavést potenciály. Pokud je vektorové pole o intenzitě E (→ r ) nevírové, platí − → rot E = 0 (68) − → Podle rovnice (67) existuje pro toto pole skalární funkce φ( r ) taková, že platí → − E = − grad φ

(69)

Tuto funkci nazveme (skalárním) potenciálem24 . Nevírové vektorové pole tedy můžeme popsat intenzitou nebo potenciálem. Popis pomocí potenciálu má řadu výhod, mj. je matematicky jednodušší. Potenciál φ však nemá sám o sobě fyzikální význam. Je to jen pomocná matematická veličina. To plyne ze skutečnosti, že je určen až na konstantu. V elektrostatice lze však z definovaného potenciálu φ snadno vypočíst potenciální energii náboje q v elektrostatickém poli. Hodnotu potenciální energie získáme když známý potenciál vynásobíme nábojem, W = qφ. Energie přitom patří mezi klíčové fyzikální veličiny. Neurčitost v konstantě se pak v případě nutnosti obchází volbou nulové hladiny potenciální energie.

PŘÍKLAD 13 Jednoduchým příkladem pro potenciál může být homogenní elektrické pole, které má pro jed− → noduchost směr osy Z, tj. E o = (0, 0, Eo ). Pak jsou jeho potenciál φ a potenciální energie We bodového náboje o velikosti q podle definice (69) a předchozího textu po řadě dány vztahy φ = −Eo z,

We = −qφ = −qE0 z.

(70)

Ze vzorců (70) je jasně patrná míra nejednoznačnosti určení potenciálu a potenciální energie náboje v elektrickém poli. V uvažovaném příkladu obojí závisí na zvolené poloze počátku souřadné osy Z. Dále bude ukázáno, že nejednoznačnost v určení potenciálu nebrání výpočtu velikosti práce nutné k přesunu náboje v elektrickém poli. • Nejprve uvažujme případ se kladný bodový náboj |q| působením elektrostatické síly Fz = |q|Eo přesune podél osy Z z místa o souřadnici z1 do vzdálenějšího bodu se souřadnicí z2 , při čemž z2 > z1 . Energie tohoto náboje se přesunem změnío hodnotu ∆We = −|q|Eo (z2 − z1 ) = |q|Eo (z1 − z2 ) < 0.

(71)

Potenciální energie náboje tedy poklesne. • Ve druhém případě mějme záporný náboj −|q|. Ten se působením elektrostatické síly přesune ze vzdálenějšího bodu se souřadnicí z2 , do bližšího bodu o souřadnici z1 , z1 < z2 . Jeho energie se změní o hodnotu ∆We = −(−|q|)Eo (z1 − z2 ) = |q|Eo (z1 − z2 ) < 0.

(72)

Nyní má konečná poloha souřadnici z1 a počáteční souřadnici z2 . Také v tomto případě potenciální energie náboje poklesne. Energie náboje (a tedy i elektrické soustavy tvořené nábojem a elektrostatickým polem) v obou možných případech klesne, elektrostatické pole koná práci. Zákon zachování energie není porušen, poněvadž současně musí existovat mechanická soustava, která zafixuje konečnou polohu 24

V rovnici (69) je znaménko minus. To souvisí s definicí znaménka změny energie, jak uvidíme dále.

47

náboje. Soustava může být tvořená např. myšlenou ideální pružinou, k níž je náboj připojen. Po zapnutí pole se podle znaménka náboje pružina buď roztáhne, nebo stlačí. Její potenciální energie se v každém případě zvýší. Stejně snadno jako v uvedeném příkladu skalárního potenciálu lze vysvětlit změnu potenciální energie zkušebního náboje v kulově symetrickém poli bodového náboje. Znaménko minus v definiční rovnici (69) je logicky vysvětleno. Vektorový potenciál → − → V případě, že vektorové pole B(− r ) je nezřídlové, platí → − div B = 0. − →→ Podle rovnice (66) existuje pro toto pole vektorová funkce A(− r ) taková, že platí ( ) − → − → ∂Ay ∂Ax ∂Az ∂Ay ∂Az ∂Ax B = rot A = − , − , − . ∂y ∂z ∂z ∂x ∂x ∂y

(73)

(74)

− → Veličinu (vektorovou funkci) A nazveme vektorovým potenciálem25 . Analogicky k nevírovému → − poli můžeme nezřídlové vektorové pole popsat buď vektorem B, který vystihuje jeho fyzikální → − podstatu (zejména silové působení), nebo vektorovým potenciálem A. Popis pomocí vektorového potenciálu je obvykle matematicky jednodušší. Vektorový potenciál však nemá fyzikální význam. To plyne již z faktu, že není určen jednoznačně. K vektorovému potenciálu lze přičíst gradient libovolné skalární funkce u(x, y, z) a výsledkem rotace je totéž vektorové pole. Uvedené tvrzení lze ověřit na základě následujících vzorců. ( ) − → ∂u ∂u ∇u(x, y, z) = ∂u , , ∂x ∂y ∂z ) − → ( ∂u ∂u C = Ax + ∂x , Ay + ∂u , A + z ∂y ∂z ) − → ( ∂Cz − → ∂Cy ∂Cx ∂Cx z ∂Cy B = rot C = ∂y − ∂z , ∂z − ∂C , − . ∂x ∂x ∂y Podobně jako u skalárního potenciálu lze v magnetizmu z vektorového potenciálu odvodit potenciální energii proudu v magnetickém poli.

PŘÍKLAD 14 − → Uvažujeme homogenní magnetické pole s indukcí B0 ve směru osy Z, tj. B o = (0, 0, Bo ). Vek− → torový potenciál uvažovaného pole s intenzitou B 0 je libovolné z řešení soustavy (75). ∂A

y z 0 = ∂A ∂y − ∂z ∂Az x 0 = ∂A ∂z − ∂x ∂Ay x B0 = ∂x − ∂A ∂y

(75)

K získání jednoznačného řešení soustavy není zadán dostatek informací. Chybí okrajová podmíka → − A(x0 , y0 , z0 ). 25

− → V případě magnetického pole je veličinou B magnetická indukce.

48

Se známými informacemi může být vektorový potenciál pole určen například následujícími vztahy ( ) − → 1 1 As = B0 y, − B0 x, 0 , 2 2 − → (76) A l = (B0 y, 0, 0), − → A p = (0, −B0 x, 0). − → Správnost řešení je možné ověřit dosazením vektorových potenciálů A i , i ∈ {l, p, s} do (75). Snadno se najde skalární funkce, jejíž gradient převede jeden výraz v (76) na druhý. Neurčitost ve volbě vektorového potenciálu umožňuje klást na něho další podmínku. Obvykle je to − → podmínka div A = 0. Všechny vektory ve vztahu (76) ji splňují.

4.5 Důležité vztahy vektorové analýzy V této kapitole uvedené vztahy se týkají působení operací gradient, divergence a rotace na složitější skalární či vektorové funkce. Často je to součin skalární a vektorové funkce. Relativně snadno lze některé vztahy symbolicky vyjádřit pomocí vektorového operátoru nabla, symbolicky − → onačeného jako ∇ a zavedeného v kapitole 4.2. V některých matematických a fyzikálních vztazích se vyskytuje vekorová operace − → − → − → ∂ ∂ ∂ + Ay + Az . A · grad = A · ∇ = Ax ∇x + Ay ∇y + Az ∇z = Ax ∂x ∂y ∂z

(77)

PŘÍKLAD 15 Abychom osvětlili využití zavedeného operátoru a jeho fyzikální význam, připomeňme, že vzorec − → (77) umožňuje vypočítat sílu, kterou nehomogenní magnetické pole o indukci B působí na − elementárním magnetický dipól o dipólovém momentu → m → − − → − F = (→ m · grad) B.

(78)

Uvedený vztah vyjadřuje tři rovnice pro jednotlivé složky síly: ∂Bx ∂Bx ∂Bx + my + mz , ∂x ∂y ∂z ∂By ∂By ∂By F y = mx + my + mz , ∂x ∂y ∂z ∂Bz ∂Bz ∂Bz Fz = mx + my + mz . ∂x ∂y ∂z

F x = mx

(79)

− → Významné uplatnění nachází vektorový operátor ∇ při jednoduchém odvozování důležitých − → vztahů. Aplikace vektorového operátoru ∇ musí respektovat jak pravidla pro derivování, tak pro algebraické vektorové operace. Nejjednoušší je odvození vztahu pro gradient součinu skalárních funkcí − → − → − → grad(φψ) = ∇(φψ) = ψ ∇φ + φ ∇ψ = ψ grad φ + φ grad ψ (80) Pro divergenci součinu skalární a vektorové funkce lze odvodit vztah − → → − → → − → − → − → − → − → − − div(φ A) = ∇ · (φ A) = ∇φ · A + φ ∇ · A) = A · grad φ + φ div A. 49

(81)

− → Podle pravidla o derivaci součinu dvou funkcí musíme aplikovat operátor ∇ odpovídajícím způsobem na první i na druhou funkci v součinu. Na skalární pole φ je možné aplikovat jen gradi− → ent, na vektorové pole A jen divergenci. V prvním sčítanci ve vztahu (81) je nutno respektovat, že se jedná o skalární součin dvou vektorů. Ve druhém sčítanci jde součin dvou skalárů. Připoměňme, že vztah (81) je v elektrostatice použit při odvození modelu polarizovaného dilektrika. Účinek elektrickýh dipólů lze nahradit fiktivními vázanými povrchovými a objemovými náboji. Podobným způsobem jako v předcházejících dvou případech lze odvodit vztah (82) pro rotaci součinu skalární a vektorové funkce. − → − → − → → − → − → − → − → → − − rot(φ A) = ∇ × (φ A) = ∇φ × A + φ ∇ × A) = grad φ × A + φ rot A.

(82)

Vzorec (82) se od (81) liší jen náhradou skalárního součinu vektorovým. V tomto případě je ale již nutno dbát na pořadí vektorů a s vektory prováděných operací. Vztah (82) je využit při popisu zmagnetizovaného materiálu. Magnetické dipóly lze pomocí něho nahradit vázanými povrchovými a objemovými proudy. Snadno lze za využití vztahu (15) z kapitoly 3.3 odvodit také vzorec pro dvojnásobou rotaci − → − → → − → − → − → − vektoru. Stačí v (15) provést jednoduchou substituci A = ∇, B = ∇ a C = A. Dostaneme − → → − → − → − → − → − → − − →− → → − → − →− →− → − rot(rot( A)) = ∇ × ( ∇ × A) = ∇( ∇ · A) − ( ∇ · ∇) A = ∇( ∇ · A) − ∆ A.

(83)

Odvozený vztah (83) je snad nejdůležitější v teorii elektromagnetického pole. Pomocí něho se odvozují Laplaceova, Poissonova rovnice, rovnice pro skinefekt, vlnová rovnice, rovnice pro elektromagnetické potenciály apod. Další důležitý vztah má tvar → − − → − − → → → − − → div( A × B) = B · rot A − A · rot B.

(84)

Ve dvou níže uvedených příkladech následuje demonstrace platnosti vztahů (83) a (84).

PŘÍKLAD 16 V programu Sage jsou nadefinovány funkce pro výpočet rotace vektoru (curl), divergence a Laplaceův operátor26 .

sage: var('x y z') (x, y, z) sage: a1 = function('a1',x,y,z) sage: a2 = function('a2',x,y,z) sage: a3 = function('a3',x,y,z) sage: def curl(F): assert(len(F) == 3); return vector([diff(F[2],y) - ... diff(F[1],z), diff(F[0],z) - diff(F[2],x), diff(F[1],x) - diff(F[0],y)]) sage: def divergence(F): assert(len(F) == 3); return (diff(F[0],x) + ... diff(F[1],y) + diff(F[2],z)) sage: def laplace(F): assert(len(F) == 3); return vector([diff(F[0],x,2) + ... diff(F[0],y,2) + diff(F[0],z,2), diff(F[1],x,2) + diff(F[1],y,2) + ... diff(F[1],z,2), diff(F[2],x,2) + diff(F[2],y,2) + diff(F[2],z,2)]) sage: A = vector([a1, a2, a3]) sage: L = curl(curl(A)) sage: pom(x,y,z) = divergence(A) sage: R = pom.diff() - laplace(A) 26 Symbolem …“ je níže označeno pokračování kódu na dalším řádku. Do příkazové řádky se tečky nepíší a kód ” se nerozděluje.

50

Níže jsou pro ilustraci porovnány druhé složky (index [1]) vektorů z obou stran rovnice (83). L[1] = R[1] =

∂ 2 a1 (x,y,z) ∂x∂y ∂ 2 a1 (x,y,z) ∂x∂y

− −

∂ 2 a2 (x,y,z) ∂x2 ∂ 2 a2 (x,y,z) ∂x2

− −

∂ 2 a2 (x,y,z) ∂z 2 ∂ 2 a2 (x,y,z) ∂z 2

+ +

∂ 2 a3 (x,y,z) ∂y∂z ∂ 2 a3 (x,y,z) ∂y∂z

Porovnání zbylých složek vektorů z levé a pravé strany rovnice (83) je ponecháno na čtenáři.

PŘÍKLAD 17 Za pomoci programu Sage a několika následujících algebraických úprav lze platnost vztahu (84) ověřit níže uvedeným způsobem. Při ověření jsou využity dříve nadefinované funkce pro výpočet divergence a rotace (curl) vektoru. sage: var ('x y z') (x, y, z) sage: a1 = function('a1',x,y,z) sage: a2 = function('a2',x,y,z) sage: a3 = function('a3',x,y,z) sage: b1 = function('b1',x,y,z) sage: b2 = function('b2',x,y,z) sage: b3 = function('b3',x,y,z) sage: A = vector([a1, a2, a3]) sage: B = vector([b1, b2, b3]) sage: L = divergence(A.cross_product(B)) sage: R = B.dot_product(curl(A)) - A.dot_product(curl(B))

Platnost prověřované rovnice lze provést porovnáním L a R. ∂b2 ∂b1 ∂b3 ∂b1 ∂b2 ∂a2 ∂a3 3 L = ( ∂b ∂y − ∂z ) · a1 + ( ∂z − ∂x ) · a2 − ( ∂y − ∂x ) · a3 + ( ∂z − ∂y ) · b1 − ∂a3 ∂a1 ∂a2 1 −( ∂a ∂z − ∂x ) · b2 + ( ∂y − ∂x ) · b3 ∂a1 ∂a2 ∂a2 ∂a3 ∂a3 ∂b1 1 R = b3 · ∂a ∂y − b2 · ∂z − b3 · ∂x + b1 · ∂z + b2 · ∂x − b1 · ∂y − a3 · ∂y + ∂b1 ∂b2 ∂b2 ∂b3 ∂b3 +a2 · ∂z + a3 · ∂x − a1 · ∂z − a2 · ∂x + a1 · ∂y

− → Z používaných vztahů zbývají ještě dva, které nelze ani s pomocí operátoru ∇ na papíře odvodit příliš snadno. Pro gradient skalárního součinu dvou vektorů platí vztah − → − − → − → − → − → → → − − → → − − → grad( A · B) = ( A × rot B) + ( B × rot A) + ( B · grad) A + ( A · grad) B.

(85)

− → Pomocí operátoru ∇ lze tento vztah přepsat na vztah → − − → → → → − − → − → − − → − − − →− − → → − → − →− →→ ∇( A · B) = ( A × ( ∇ × B)) + ( B × ( ∇ × A) + ( B · ∇) A + ( A · ∇) B.

(86)

Pro rotaci vektorového součinu dvou vektorů platí vztah − → − − → → − → − − → → → − − → − → − → rot( A × B) = ( A div B) − ( B div A) + ( B · grad) A − ( A · grad) B.

(87)

− → S využitím operátoru ∇ dostane uvedený vztah tvar − → → − →− − → − → → − → − →− − → − → − → → − → → →− − − ∇ × ( A × B) = ( A ∇ · B) − ( B ∇ · A) + ( B · ∇) A − ( A · ∇) B.

51

(88)

PŘÍKLAD 18 Platnost rovnic (86) a (88) jde dokázat postupem analogickým důkazu vztahu (84). Jen je − → − → →− potřeba dodefinovat a použít funkci pro vektorovou operaci ( A · ∇) B. Porovnávání je víc, protože výsledkem na obou stranách rovnic (86) a (88) jsou vektory a souhlasit musí všechny jejich složky. sage: def magie(F,G): assert(len(F) == 3); assert(len(G) == 3); return vector([F[0]*diff(G[0],x) + F[1]*diff(G[0],y) + F[2]*diff(G[0],z), F[0]*diff(G[1],x) + F[1]*diff(G[1],y) + F[2]*diff(G[1],z), F[0]*diff(G[2],x) + F[1]*diff(G[2],y) + F[2]*diff(G[2],z)])

4.5.1 Shrnutí Na závěr této části shrneme nejdůležitější uvedené vztahy jednak v plné formě, která využívá standardní symboly pro vektorové operace, jednak ve zkrácené formě s vektorovým operátorem → − ∇. V plné formě mají vzorce tvar: − → − → → − rot(rot( A) = grad(div A) − ∆ A, − → → − − → div(φ A) = A · grad φ + φ div A, − → − → − → rot(φ A) = grad φ × A + φ rot A, → − − → − − → → → − − → div( A × B) = B · rot A − A · rot B, grad(φψ) = ψ grad φ + φ grad ψ, − → − − → − → − → − → → → − − → → − → − grad( A · B) = ( A × rot B) + ( B × rot A) + ( B · grad) A + ( A · grad) B, − → − − → → − → − − → → → − − → − → − → rot( A × B) = ( A div B) − ( B div A) + ( B · grad) A − ( A · grad) B. − → Využitím vektorového operátoru ∇ přejdou vztahy do podoby rovnic: → − → − → → − − → − − → → − → − →− ∇ × ( ∇ × A) = ∇( ∇ · A) − ( ∇ · ∇) A, → − → → − → − → − → − − ∇ · (φ A) = ∇φ · A + φ ∇ · A), − → − → − → − → − − → → ∇ × (φ A) = ∇φ × A + φ ∇ × A), → − → − → − → − → − → − → − → − → − ∇ · ( A × B) = B · ∇ × A − A · ∇ × B, → − → − − → ∇(φψ) = ψ ∇φ + φ ∇ψ, → − − → → → → − − → − → − − → − − − →− − → → − → − →− →→ ∇( A · B) = ( A × ( ∇ × B)) + ( B × ( ∇ × A) + ( B · ∇) A + ( A · ∇) B, − → → − →− − → − → → − → − →− − → − → − → → − → → →− − − ∇ × ( A × B) = ( A ∇ · B) − ( B ∇ · A) + ( B · ∇) A − ( A · ∇) B.

4.6 Matematický popis vektorového pole Níže v této kapitole naleznete matematickou formulaci a několik ekvivalentních vyjádření popisu nevírových, popřípadě nezřídlových vektorových polí. Předtím však uvedeme, některé méně známé matamatické symboly, které se zde používaji. • Symbolem Rn je označen prostor reálných čísel dimenze n. V našem případě uvažujeme prostor dvourozměný (zřídka) nebo trojrozměrný. Proto n = 2, 3. • Oblast v prostoru Rn označíme Ω, tedy Ω ⊂ Rn (Ω je podmnožina Rn ). Tato oblast může být i degenerovaná, o nižším počtu rozměrů, např. Ω = R1 . • Skalární funkci f popíšeme tímto způsobem: f : Ω → R1 , Ω ⊂ Rn . Předpis říká, že definičním oborem fukce f je oblast Ω a funkce f ji zobrazuje do prostoru R1 . Je to prostor dimenze 1, poněvadž jde o skalární funkci. Výsledkem funkce f je skalár, tj. jedno reálné číslo. 52

• Pro vektorovou funkci A můžeme psát analogicky: A : Ω → Rn . V našem případě jde o zobrazení do prostoru dimenze 2, nebo 3, podle volby. • Symbol ∃φ znamená, že existuje v tomto případě funkce φ s vlastnostmi, které následují dále za tímto výrazem. • Symbol ∀c říká, že všechny objekty c, např. křivky, mají vlastnost, která následuje.

PŘÍKLAD 19 Jako příklad použití uvedených symbolů uvedeme definici spojitosti reálné funkce jedné reálné proměnné. Slovně je např. takto: Reálná funkce f jedné reálné proměnné x, tj. f (x), je spojitá v bodě x0 , jestliže ke každému ϵ > 0 existuje δ > 0 tak, že pro všechna |x − xo | < δ platí |f (x) − f (x0 )| < ϵ. Pomocí výše uvedených symbolů lze definici spojitosti psát např. takto: f : Ω → R1 , ∀ϵ > 0

Ω ⊂ R1 je spojitá v bodě x0 , pokud

∃δ > 0

∀(|x − xo | < δ) :

|f (x) − f (x0 )| < ϵ.

Zápis je stručný a přehledný. Zejména pro začátečníka může však být matoucí.

Potenciální (konzervativní, nevírové) pole − → Vektorové pole A : Ω → Rn , kde n = 2, 3 a Ω ⊂ Rn , se nazývá potenciální na Ω, jestliže − → integrace A po křivkách ležících v Ω nezávisí na cestě, tj. pokud pro křivky c1 a c2 , ležící v Ω a mající stejný počáteční a stejný koncový bod, platí ∫ ∫ − → − − → − → → Ad l = Ad l . (89) (c1 )

(c2 )

Věta 1.: − → Pokud funkce A : Ω → Rn , kde n = 2, 3 a Ω ⊂ Rn , kterou je zadáno vektorové pole, má spojité parciální derivace na Ω, následující tvrzení jsou ekvivalentní: → − 1. Vektorové pole A je na Ω potenciální, → − → → − − . 2. ∇ × A = rot A = 0 v Ω, ∫− → − → 3. Ad l = 0, ∀c jednoduše uzavřené hladké křivky v Ω, c

→ − − → 4. ∃φ : Ω → R taková, že ∇φ = gradφ = A na Ω. Příkladem konzervativního vektorového pole je například gravitační pole.

53

Nezřídlové (solenoidální) pole − → Vektorové pole A : Ω → R3 na oblasti Ω ⊂ R3 se nazývá nezřídlové, jestliže pro libovolné dvě hladké, uzavřené plochy S1 , S2 , které leží v Ω, mají stejný okraj (ohraničující křivku (C), viz obr. 20) a jsou souhlasně orientovány vzhledem k tomuto okraji, platí ∫ ∫ − → − − → → → − Ad S = Ad S . (90) (S1 )

(S2 )

Obrázek 20: Ilustrace k definici solenoidálního pole. Věta 2.: − → Pokud funkce A : Ω → R3 , a Ω ⊂ R3 , kterou je zadáno vektorové pole, má spojité parciální derivace na oblasti Ω, následující tvrzení jsou ekvivalentní: → − 1. Vektorové pole A je na Ω nezřídlové, → → − → − − . 2. ∇ · A = div A = 0 na Ω, ∫− → − → 3. Ad S = 0 pro libovolnou uzavřenou plochu v S ⊂ Ω, S

→ − → − → − → − 4. ∃ B : Ω → R3 taková, že A = ∇ × B = rot B. − → Vektorové pole B s vlastností 4 se nazývá vektorový potenciál“. ”

54

5 Základy tenzorového počtu Tato část logicky navazuje na dvě předchozí, pojednávající o vektorech. Pomocí tenzorů druhého řádu lze popsat obecnou závislost mezi dvěma vektory. Tenzory vyšších řádů pak popisují ještě složitější závislosti. Zde se soustředíme pouze na základy tenzorového počtu, a to především na tenzory druhého řádu. Ty popisují např. elektrické a magnetické vlastnosti krystalů, transportní jevy v anizotropním prostředí apod. Jsou základem teorie pružnosti, protože deformace a elastické napětí jsou tenzory druhého řádu. Tenzorový počet je rozsáhlá disciplína. Je popsána v řadě monografií, za zdařilou považujeme například učebnici [15]. V této části knihy budou poskytunuty pouze základní informace. Důraz bude kladen na pochopení tenzorů v elasticitě, která začíná hrát důležitou roli v moderních aplikacích včetně mechatronických.

5.1 Tenzor elektrické vodivosti Výklad tenzorů začneme materiálovým tenzorem elektrické vodivosti γ, na kterém jsou tenzorové vlastnosti vidět nejlépe. Tenzor elektrické vodivosti popisuje vztah (91) mezi vektorem intenzity − → − → elektrického pole E a vektorem proudové hustoty i . − → − → i = γE

(91)

− → − → V izotropním prostředí mají vektor proudové hustoty i a intenzity elektrického pole E stejný směr, viz obr. 21(a). Zásadní odlišnost anizotropního prostředí (vodivého krystalu) je − → − → v tom, že směry vektorů i a E nejsou shodné, viz obr. 21(b). Navíc se úhel mezi vektory mění, − → pokud se mění směr budící složky, tedy vektoru intenzity elektrického pole E .

(a)

(b)

Obrázek 21: Vztah mezi proudovou hustotu a elekrickým polem v izotropním (21(a)) a anizotropním (21(b)) protředí. Podle obr. 21(b), může v anizotropním prostředí elektrické pole ve směru souřadné osy obecně vyvolat všechny tři složky proudové hustoty. V takovém případě musí Ohmův zákon (91) tuto skutečnost respektovat. Proto má pro anizotropní prostředí složitější tvar       E1 γ11 γ12 γ13 i1  i2  =  γ21 γ22 γ23  ·  E2  , (92) E3 γ31 γ32 γ33 i3 kde fyzikální veličina

 γ11 γ12 γ13 γ =  γ21 γ22 γ23  γ31 γ32 γ33 

(93)

je tenzor elektrické vodivosti. Je to tenzor druhého řádu, má devět složek. Pomocí energetických úvah lze dokázat, že je symetrický γij = γji , takže nezávislých složek je nejvýše šest. 55

Pomocí indexů se Ohmův zákon v diferenciálním tvaru (92) pro anizotropní prostředí dá napsat jako 3 ∑ ii = γij Ej . (94) j=1

Předchozí zápis se zjednodušuje pomocí Einsteinova sumačního pravidla. Přes dvakrát se opakující index se sčítá od jedné do tří. (95)

ii = γij Ej Jednoduchou praktickou aplikaci ukážeme na příkladu.

PŘÍKLAD 20 − → − → Nechť elektrické pole o intenzitě E má směr osy X1 , tj. E = (E1 , E2 , E3 ) = (E, 0, 0). Pak z kterékoliv formulace Ohmova zákona v tenzorovém tvaru ((92), (94), (95)) zjitíme, že toto elekrické pole vyvolá všecny tři složky proudové hustoty i1 = γ11 E1 = γ11 E,

i2 = γ21 E1 = γ21 E,

i3 = γ31 E1 = γ31 E

(96)

Proud má tedy obecně nenulovou složku ve směru každé ze tří souřadných os. Teoretický ověřovací experiment je zachycen na následujícím obrázku.

Obrázek 22: Elektrické schéma experimentu. V případě, že je tyčinka z vodivého anizotropního krystalu, ideálním ampérmetrem A1 změříme velký hlavní proud I1 , avšak ideální ampérmetr A2 ukáže obecně malý, ale nenulový proud, I2 . Vodivost γ izotropního protředí je jako skalár speciálním případem tenzoru druhého řádu. Použitím koeficientu Kroneckerovo-delta, který je definován vztahy27 δij = 1, pro i = j, δij = 0, pro i ̸= j,

(97)

můžeme Ohmův zákon pro izotropní prostředí psát ve tvaru ii = γδij Ej .

(98)

27 V tomto vztahu nemůžeme psát δii = 1, poněvadž podle Einsteinova sumačního pravidla, které zde požíváme, se podle dvakrát opakujícího se indexu i sčítá od jedné do tří, tedy δii = δ11 + δ22 + δ33 = 1 + 1 + 1 = 3

56

Tenzor vodivosti izotropního prostředí má tedy tvar   γ 0 0 γ (i) =  0 γ 0  = γδij . 0 0 γ

(99)

5.2 Tenzory v mechanice pružných těles Uvažujme pružné těleso vystavené působení povrchových sil. Vlivem silového působení může dojít ke změně tvaru pružného tělesa (deformace). Malé změny tvaru pružných těles jsou obvykle vratné (elastická deformace). V takovém případě se těleso po konci působení sil navrátí do původního tvaru. Působení plošné síly se při elastické deformaci projevuje vznikem napětí (síla vztažená na plochu) uvnitř tělesa. Tato kapitola přináší ukázku tenzorového popisu napětí a napětím vyvolané deformace tělesa. Vztah mezi tenzorem deformace a tenzorem napětí vyjadřuje takzvaný Hookeův zákon nastíněný v kapitole 5.4. Kromě tenzorů deformace a napětí je dále v kapitole představen tenzor otočení tělesa. Kombinace tenzoru otočení a tenzoru deformace umožňuje popsat prostý střih. 5.2.1 Tenzor deformace Tenzor deformace je tenzorem druhého řádu. Neprovedeme systematický výklad, jen znímíme, → že při silovém působení na poddajné těleso se bod o polohovém vektoru − r přesune do místa − → → − ′ o polohovém vektoru r . Pracuje se ale s vektorem posunutí u , což je rozdíl obou vektorů, tedy − → → → r′=− r +− u, ri′ = ri + ui , i = 1, 2, 3, (100) − → − → u =Sr. Pro malá posunutí ui se pak projeví lineární deformace (viz např. [15]), která je popsána tenzorem deformace ( ) ∂ui 1 ∂uj + . (101) Sij = 2 ∂xi ∂xj Z definice je zřejmé, že tenzor deformace je symetrický, Sij = Sji . Má tedy maximálně šest nezávislých složek. Stejně jako u tenzoru napětí mají diagonální a nediagonální složky odlišný význam. Význam diagonálních složek je na obr. 23. Krátká tyč malého průřezu původní délky ∆x1 ve směru osy X1 se prodlouží o délku ∆u1 ve směru téže osy. Podle definice je její deformace ( ) 1 ∆u1 ∆u1 ∆u1 S11 = + = (102) 2 ∆x1 ∆x1 ∆x1 Limitním přechodem ∆x1 → 0 získáme diagonální složku S11 tenzoru deformace (101). S11 =

∂u1 ∂x1

(103)

Analogicky to platí pro ostatní diagonální složky, představují tedy relativní prodloužení ve směru os souřadné soustavy. K vysvětlení nediagonálních složek tenzoru deformace (101) použijeme obr. 24, který popisuje deformaci v rovině X1 X2 . Malý čtverec se stranami ∆x1 a ∆x2 (∆x1 = ∆x2 ) se působením blíže nespecifikované síly deformuje v kosočtverec tak, že vrchol původně na ose X1 se posune o vzdálenost ∆u2 ve směru osy X2 a vrchol původně na ose X2 se posune o vzdálensot ∆u1 ve směru osy X1 . Hrany kosočtverce svírají s původními hranami čtverce malý úhel δ/2. Poněvadž tento úhel je malý, platí pro něho vztahy δ . δ . ∆u2 ∆u1 = tg = = . 2 2 ∆x1 ∆x2 57

(104)

Obrázek 23: Podélná deformace

Obrázek 24: Střižná (příčná, smyková) deformace. Z těchto vztahů plyne 1 2

(

∆u2 ∆u1 + ∆x1 ∆x2

)

. δ = . 2

Limitním přechodem ∆x1 → 0 a ∆x2 → 0 při zvážení vztahu (105) dostaneme vztah ( ) δ 1 ∂u2 ∂u1 + = . S12 = 2 ∂x1 ∂x2 2

(105)

(106)

Přibližný vztah (105) přejde v limitě na přesný vztah (106). Odtud plyne fyzikální význam nediagonální složky S12 tenzoru deformace. Složka S12 udává polovinu úhlu28 , o který se změní původní pravý úhel mezi osami X1 a X2 . Osy, o které se jedná, jsou obsaženy v indexech složky S12 . Pokud je složka S12 kladná, pravý úhel se zmenší na π/2 − δ = π/2 − 2 · S12 , viz obr. 24. Pro zápornou složku se pravý úhel zvětší na π/2 + |δ| = π/2 + 2 · |S12 |. Tuto deformaci nazýváme střižnou nebo smykovou. Analogický význam mají dvě zbývající nezávislé nediagonální složky, např. S23 je střih v rovině X2 X3 , udává tedy, oč se změní pravý úhel v této rovině. Složka S13 pak popisuje střih v rovině X1 X3 . Obecně složka Sij , i ̸= j, udává střih v rovině Xi Xj . 5.2.2 Tenzor napětí Odvození a vysvělení uváděných vztahů lze najít v učebnicích, zde se soustředíme na fyzikální význam a aplikace. Zjednodušeně řečeno, elastické napětí se projeví silou působící na malou plochu. K jeho určení je tedy nutno uvést složky síly a popsat orientaci (malé) plochy. Orientace plochy se nejjednodušeji popíše jednotkovým vektorem normály k ní29 . Elastické napětí je faktor v linárním − → → − − n plochy ∆ S vztahu (107) mezi dvěma vektory, jeden určuje složky síly ∆ F a druhý normálu → na níž síla působí. 28

Úhly udáváme v radiánech. Na stupně je nutno je přepočíst. Plocha má dvě strany, proto vektor normály má určený jen směr, nikoliv smysl. Tato nejednoznačnost v praktických aplikacích nevadí, poněvadž se vyšetřuje působení na objemové elementy a pak je smysl normály jednoznačný, ven z elementu. 29

58

Známe-li složky tenzoru napětí Tij , můžeme spočítat sílu působící na malou rovinnou plochu − → − ∆S = |∆ S | s vektorem normály → n = (n1 , n2 , n3 ) podle vztahu ∆Fi = Tij nj ∆S.

(107)

Obr. 25 to vysvětluje názorně.

Obrázek 25: Tenzor napěti a síla Výklad vlastností tenzoru napětí provedeme nejprve za těchto předpokladů: 1. Napětí je homogenní ve vyšetřované oblasti. 2. Vyšetřovaná oblast je ve statické rovnováze. 3. Neexistují objemové síly a momenty. Tenzor napětí



 T11 T12 T13 T =  T21 T22 T23  T31 T32 T33

(108)

má tedy složky o dvou indexech. První index složky tenzoru napětí pak udává směr působící síly, index 1 sílu působící ve kladném směru osy X1 , atd. Druhý index popisuje normálu plochy, na kterou bude síla působit. Např. druhý index o hodnotě 1 znamená, že se jedná o plochu kolmou k ose X1 , jejíž normála směřuje v kladném směru osy X1 , složky napětí o druhém indexu rovném 2 se vztahují k ploše kolmé k ose X2 s kladnou normálou v jejím směru, atd. Toto platí zcela obecně pro diagonální i nediagonální složky tenzoru napětí. Z orientace síly a normály vyplývá též smysl kladných a záporných hodnot složek tenzoru napětí. Složka napětí je kladná, pokud 1. normála plochy směřuje v kladném směru příslušné osy a rovněž síla je v kladném směru odpovídající osy, nebo 2. normála plochy směřuje proti směru příslušné osy a rovněž síla je proti směru odpovídající osy. Ve zbývajícíh dvou případech je znaménko záporné: 1. Normála plochy je v záporném směru příslušné osy a síla je v kladném směru odpovídající osy, nebo 2. normála plochy směřuje proti směru příslušné osy, zatímco síla je ve směru odpovídající osy.

59

Podle tohoto popisu se odlišně chovají diagonální a nediagonální složky tenzoru napětí, a proto je probereme odděleně. Význam diagonálních složek by měl vyplynout z obr. 26. Na krychli se stěnami kolmými k souřadným osám působí napětí kolmá k těmto plochám. Nazýváme je normálová napětí. Normály k plochám směřují ven z krychle. Na protilehlé strany krychle musí působit tatáž napětí. Jinak by nebyla dodržena statická rovnováha sil a krychle by se vlivem výsledné síly pohybovala rovnoměrně zrychleně. Podle obr. 26 na plochu s normálou ve směru osy X1 působí napětí T11 ve směru téže osy. Je kladné, poněvadž jak normála, tak síla, mají kladný směr osy X1 . Na protilehlou plochu působí síla ve směru záporné osy X1 , normála k ploše je rovněž ve směru záporné osy X1 , normálové napětí T11 je tedy kladné i na této ploše.

Obrázek 26: Význam diagonálních, složek tenzoru napětí, normálová napětí. Obdobně lze vysvětlit normálovou složku napětí T22 . Na stěnu krychle s vnější normálou ve směru kladné osy X2 působí napětí ve směru téže osy. Napětí na opačné straně má rovněž kladné znaménko, poněvadž jak normála ke stěně, tak síla na ni působící, směřují podél záporné osy X2 . Analogicky lze popsat normálovou složku napětí T33 . Kladné napětí, jak je tomu na obr. 26 působí ven z krychle. Snaží se tedy krychli roztáhnout. Záporná hodnota napětí znamená, že síla působí směrem dovnitř krychle, snaží se tedy krychli stlačit, zmenšit její rozměry. Toto je běžně užívaná konvence, ve speciálních případech tomu však může být obráceně. K vysvětlení nediagonáloních složek tenzoru napětí slouží obr. 27. Zde je pro přehlednost jen rovina X1 X2 souřadné soustavy. Na plochu s vnější normálou ve směru kladné osy X1 působí síla ve směru osy X2 . Leží v této ploše, její účinek se projeví jako tečné napětí T21 . Je kladné, poněvadž normála i síla jsou ve směru kladných os. Protilehlá strana má zápornou normálu, takže tečná síla musí působit ve směru záporné osy X2 , aby složka tečného napětí T21 zůstala kladná. Složka tečného napětí T21 vytváří silovou dvojici, viz obr. 27, která se snaží krychlí otáčet. Kromě této složky je zde složka T12 . Zde působí síla na stěny kolmé k ose X2 . Na horní ploše má tato tečná síla směr kladné osy X1 , aby složka T12 tenzoru napětí byla kladná. Na dolní ploše se zápornou normálou musí pak tečná síla směřovat proti kladné ose X1 . Podle obr. 27 složka T12 tenzoru napětí vytváří silovou dvojici, která se snaží krychli otáčet v opačném smyslu než složka T21 . Aby krychle zůstala v rovnováze, musí být tyto dvě složky stejné, T12 = T21 . Poněvadž podobnou úvahu lze udělat i pro zbývající dvě dvojice nediagonálních složek, docházíme k závěru, že tenzor napětí (108) je za výše uvedených zjednodušujících

60

Obrázek 27: Význam nediagonálních složek tenzoru napětí, tečná napětí. předpokladů symetrický, Tij = Tji ,

i ̸= j.

(109)

Má tedy nanejvýš šest nezávislých složek, v nejobecnějším případě existují tři nenulová normálová a tři nenulová tečná napětí.

PŘÍKLAD 21 Nyní uvažujeme případ, že nejsou splněny všechny předpoklady uvedené výše, tj. elastické napětí se mění od bodu k bodu a působí objemové síly a momenty. Těleso je ale stále ve stavu statické rovnováhy. Jako objemovou sílu budeme předpokládat tíhu. Pak na elementární objem dV působí síla dFi = ρgi dV . Hustota této síly je pak dána vztahem fi = ρgi ,

(110)

→ kde ρ je hustota látky a − g = (g1 , g2 , g3 ) je vektor tíhového zrychlení. Z podmínky pro rovnováhu sil ve směru souřadných os lze odvodit obecnou podmínku rovnováhy, viz např. [15]: ∂Tij + fi = 0. ∂xj

(111)

Z podmínky pro rovnováhu momentů, za předpokladu, že neexistují objemové momenty, se odvodí symetričnost tenzoru napětí (109). V mechanice objemové momenty neexisují, proto se podmínky (109) chápou jako obecné. Na závěr pojednání o tenzoru napětí uvedeme tři jednoduché praktické případy, jednosměrnou napjatost σ ve směru osy X3 , tečné napětí τ v rovině kolmé k ose X2 a vnější tlak p. Poněvadž vnější tlak působí směrem do tělesa, mají jeho hodnoty znaménko minus.       −p 0 0 0 0 τ 0 0 0 (112) T (σ) =  0 0 0  , T (τ ) =  0 0 0  , T (p) =  0 −p 0  . 0 0 −p τ 0 0 0 0 σ Tlak můžeme vyjádřit i pomocí symbolu Kroneckerovo delta (97) Tij = −pδij .

61

(113)

5.2.3 Tenzor pootočení Vedle deformace a napětí lze tenzorem druhého řádu popsat pootočení,   ω11 ω12 ω13 ω =  ω21 ω22 ω23  , ω31 ω32 ω33 Složky tenzoru pootočení jsou definovány vztahem ( ) 1 ∂uj ∂ui ωij = − . 2 ∂xi ∂xj

(114)

(115)

Na rozdíl od tenzoru deformace zde záleží v případě tenzoru pootočení na pořadí indexů. První index i u složky ωij se v první parciální derivaci vztahuje k souřadnici, xi , i ∈ {1, 2, 3}, druhý index j přináleží ke složce posunutí, uj , j ∈ {1, 2, 3}. U druhé parciální derivace ve vzorci (115) jsou indexy prohozené. Z těchto vztahů plyne, že tenzor pootočení je antisymetrický ωij = −ωji .

(116)

Diagonální složky tenzoru pootočení jsou nulové, obsahuje tedy jen tři nezávislé složky. Význam složek tenzoru pootočení vyplývá z obr. 28. Analogicky jako u střižné deformace se vychází ze čtverce se stranami rovnoběžnými s osami souřadné soustavy. Strany mají délky ∆x1 a ∆x2 (∆x1 = ∆x2 ). Čtverec se pootočí tak, že vrchol původně na ose X1 se posune o vzdálenost ∆u2 ve směru osy X2 a vrchol původně na ose X2 se posune o vzdálenost ∆u1 proti směru osy X1 , tedy ∆u1 = −∆u2 . To je zásadní rozdíl proti střižné deformaci.

Obrázek 28: Otočení objektu kolem osy X3 . Pro malé pootočení γ platí ∆u1 . . ∆u2 γ = tg γ = =− . ∆x1 ∆x2

(117)

U posledního výrazu (∆u1 /∆x2 ) musí být znaménko minus, poněvadž kladnému úhlu γ odpovídá záporné posunutí ∆u1 . Ze vztahů (117) snadno odvodíme ( ) 1 ∆u2 ∆u1 . − = γ. (118) 2 ∆x1 ∆x2 Analogicky jako u střižné defromace limitním přechodem ∆x1 → 0 a ∆x2 → 0 při zvážení vztahu (118) dostaneme vztah ( ) 1 ∂u2 ∂u1 ω12 = − . (119) 2 ∂x1 ∂x2 62

Odtud plyne fyzikální význam složky ω12 tenzoru otočení. Udává úhel, o který se objekt v rovině X1 X2 pootočí kolem osy X3 . Dochází k pootočení kolem té osy, jejíž index u složky tenzoru chybí. Smysl pootočení určují indexy, v tomto případě to je od osy X1 k ose X2 , pokud je složka ω12 kladná. Připomeňme, že kladná hodnota je pro pootočení proti smyslu pohybu hodinových ručiček. Analogický význam má složka ω23 , udává pootočení kolem osy X1 a to od osy X2 k ose X3 , pokud je složka ω23 kladná. Úhel pootočení kolem osy X2 udává složka ω31 a nikoliv složka ω13 . O tom se snadno přesvědčíme z obrázku 29, ve kterém je zachována pravotočivá soustava30 . U levotočivé soustavy změní složka ω31 znaménko. Kladná hodnota ω31 u pravotočivé soustavy znamená pootočení od osy X3 k ose X1 . Hodnota ω13 = −ω31 pak musí být záporná. Jedná se o pootočení od osy X1 k ose X3 (viz obr. 29).

Obrázek 29: Otočení objektu kolem osy X2 . Na základě předchozího rozboru má tenzor pootočení v obecném případě tvar   0 ω12 −ω13 0 ω23  . ω =  −ω21 ω31 −ω32 0

(120)

Má tedy tři nezávislé složky, tedy stejný počet jako má vektor. Obvykle se také pootočení popisuje jako vektor. V našem případě by měl složky31 − → ω = (ω1 , ω2 , ω3 ) = (ω23 , ω31 , ω12 ).

(121)

Složky tenzoru pootočení jsou ve vektorové reprezentaci (121) v cyklickém pořadí, což umožní snadné zapamatování (např. ω2 = ω31 , pořadí indexů je cyklické 2, 3, 1). Zápis (121) však platí jen pro pravotočivou soustavu, jak jsme již na příkladech ukázali. Na závěr této části připomeneme, že kombinace deformace a pootočení je v technické praxi běžná. Na obr. 30 je zakreslen prostý střih. Čtverec se v případě prostého střihu deformuje tím způsobem, že jeho strana přilehlá k ose X1 zůstává pevná a strana orientovaná původně ve směru osy X2 se pootočí o úhel −δ. 30

Osa X2 směřuje vzhůru. Používáme stejný symbol ω pro tenzor i vektor pootočení i jejich složky. Nemělo by však docházet k nedorozumění. 31

63

Obrázek 30: Prostý střih jako kombinace deformace a pootočení. Platí jenom v případě malých deformací. Prostý střih můžeme popsat nesymetrickým tenzorem32   0 0 0  δ 0 0 , 0 0 0 kde δ =

∂u1 ∂x2 .

(122)

Tento tenzor rozložíme na symetrickou a antisymetrickou část, 

  0 0 0 0  δ 0 0 = δ 2 0 0 0 0

  0 0 − 2δ 0 0  +  2δ 0 0 0 0 0 δ 2

 0 0 . 0

(123)

Z porovnání se vztahem (106) vidíme, že symetrická část představuje střižnou deformaci o úhel S12 = S21 = δ/2. Pro antisymetrickou část plyne ze vztahu (120), že jde o pootočení o úhel ω12 = −ω21 = −δ/2. Obr. 30 situaci názorně objasňuje. Nejprve se strany čtverce odchýlí v důsledku deformace od souřadných os o úhel δ/2 a −δ/2. Původně pravý úhel mezi stranami čtverce se zmenší o δ. Pak následuje pootočení deformovaného čtverce jako celku kolem osy X3 o úhel ω12 = −ω21 = −δ/2, tedy ve směru pohybu hodinových ručiček. Pootočení je od osy X2 k ose X1 , poněvadž úhel ω21 = δ/2 je kladný.

5.3 Vektory radiální a axiální V předchozí části jsme definovali vektor pootočení a upozornili, že tato definice platí jen v pravotočivé souřadné soustavě. Při přechodu k levotočivé soustavě mění jeho složky znaménko, což jsme ukázali pro pootočení kolem souřadných os. Vektor pootočení či rotace obecně řadíme mezi axiální vektory. Je kolmý k rovině otáčení a je osou otáčení. Otáčení má vůči němu axiální (osovou) symetrii. Mezi axiální vektory patří vektorový součin, který je definovan pro pravotočivou soustavu. Jak se podrobně dokazuje v učebnicích, např. [15], vektorový součin je vlastně antisymetrický tenzor druhého řádu, jehož složky při přechodu k levotočivé soustavě mění znaménko. Axiálním vektorem je tedy např. mechanický moment nebo magnetizace. Pokud se vektor chová stejně v obou soustavách (pravo- i levo-točivé), nazýváme jej radiálním. Jimi jsou (až na výjimky) všechny ostatní vektory. Z uvedeného vyplývá nutnost správně popsat souřadnou soustavu. Při chybném uvedení jejích os nemusí vztahy pro vektorový součin platit. 32

Není to ani tenzor deformace, ani tenzor pootočení

64

5.4 Tenzory vyššího řádu Tyto tenzory jsou již jen materiálové. Důležitým tenzorem třetího řádu je například piezoelektrický tenzor, který popisuje piezoelektrický jev. Dochází k němu zejména u některých krystalů. U přímého piezoelektrického jevu se mechanickým namáháním krystal elektricky polarizuje, což se makroskopicky projeví tím, že mezi vhodně umístěnými elektrodami naměříme elektrické napětí. U nepřímého piezoelektrického jevu dochází působením elektrického pole k deformaci krystalu. Piezoelektrický jev nastává mezi elektrickou veličinou (intenzita elektrického pole nebo elektrická indukce) popsanou vektorem a elastickou veličinou (elastické napětí nebo deformace), které jsou popsány tenzorem druhého řádu. K popisu piezoelektrického jevu je tedy nutný tenzor třetího řádu. Pro ilustraci uvedeme nepřímý piezoelektrický jev, kdy se působením elektrického pole budí deformace. Příslušné materiálové rovnice mají tvar Sij = dkij Ek ,

(124)

kde Ek jsou složky vektoru intenzity elektrického pole a Sij složky tenzoru deformace. Složky piezoelektrického tenzoru třetího řádu jsou dkij . Celkový počet je 27, nezávislých je maximálně 18. Tenzory čtvrtého řádu jsou v technické praxi velmi běžné. Popisují například vztah mezi elastickým napětím a deformací, tedy zobecněný Hookův zákon. Jak napětí, tak deformace, jsou tenzory druhého řádu, vztah mezi nimi musí být tenzorem čtvrtého řádu. Zobecněný Hookův zákon má dvojí tvar podle volby nezávisle proměnné. Působí-li napětí Tkl , pak pro deformaci Sij platí Sij = sijkl Tkl , (125) kde sijkl jsou elastické koeficienty. Je-li zadána deformace Skl , pak potřebné napětí Tij určíme ze vztahů Tij = cijkl Skl , (126) kde cijkl jsou elastické moduly. Oba elastické tenzory mají 81 složek, nezávislých je však maximálně 21. Jsou-li známy moduly, určíme z nich koeficienty a naopak.

65

66

6 Analogie Naprosto rozdílné fyzikální jevy mohou být popsány formálně stejnými matematickými rovnicemi a tedy řešeny stejnými metodami. Odlišné fyzikální veličiny, které v těchto rovnicích vystupují na stejných místech, nazýváme analogické a mluvíme o analogii jevů popsaných těmito rovnicemi. Analogií je ve fyzice celá řada. Zde se zmíníme o dvou z nejvýznamnějších. Jednak je to analogie mezi mechanickými (případně akustickými) a elektrickými veličinami, jednak se jedná o analogii jevů popsaných rovnicemi s potenciálovým polem. V prvním případě mluvíme o elektromechanické či elektroakustické analogii. Ve druhém případě se jedná o analogii rovnic popisujících transportní jevy.

6.1 Elektromechanická analogie Tato analogie je velmi rozsáhlá a široce využívaná, jak dosvědčují monografie [12], či novější kniha [22], která je však této analogii věnována pouze z části. Zde uvedeme jen základní výsledky, abychom získali představu o jejím využití. Elektromechanická analogie je založena na tom, že v zákonech a rovnicích pro mechanické a elektrické soustavy si navzájem odpovídají tyto veličiny • mechanická síla F a elektrické napětí U , • mechanická rychlost v a elektrický proud I, • mechanické posunutí x a elektrický náboj q. Jako příklad uvedeme rovnici pro elektrický seriový kmitavý obvod tvořený cívkou o indukčnosti L, odporem velikosti R a kondenzátorem o kapacitě C.

PŘÍKLAD 22 V RLC obvodu jsou kmity buzeny vnějším zdrojem harmonického napětí uv (t) = Uv sin(ωt + φ), kde Uv je amplituda budícího napětí, ω = 2πf je jeho úhlová frekvence (f je frekvence) a φ je fázová konstanta, která je obvykle rovna nule. Pak má difrenciální rovnice pro náboj na kondenzátoru tvar d2 q 1 dq L 2 + R + q = Uv sin(ωt + φ) (127) dt dt C

(a) RLC obvod

(b) Mechanická analogie RLC obvodu

Obrázek 31: Příklad elektromechanické analogie.

67

Nechť je mechanický oscilátor tvořen tuhým tělesem o hmotnosti m, pružinou o poddajnosti k (převrácená hodnota tuhosti) a celkovým tlumením, které formálně popíšeme odporem r. Viz druhý z výše vložených obrázků. Pro pružinu a tlumení platí vztahy x = −kFe ,

Fr = rv = r

dx , dt

(128)

kde x je výchylka, v je rychlost, Fe je elastická síla potřebná k zajištění výchylky x a Fr je síla způsobená tlumením33 . Dále předpokládáme, že vynucené mechanické kmity vyvolává vnější síla fv (t) = Fv sin(ωt + φ), kde Fv je její amplituda a význam ostatních symbolů je stejný jako u budícího elektrikého napětí. Analogie vyžaduje, aby obvod byl duální k elekrickému, tedy všechny mechanické prvky jsou spojeny paralelně. Pak lze pro tento obvod zapsat pohybovou rovnici pro posunutí d2 x dx 1 +r + x = Fv sin(ωt + φ). (129) dt2 dt k Z porovnání rovnic (127) a (129) je zřejmé, že všechny výše uvedené analogické veličiny jsou na stejným místech. Navíc zjišťujeme, že hmotnosti přísluší indukčnost, mechanickému tlumení elektrický odpor a mechanické poddajnosti elektrická kapacita. Na základě této analogie můžeme každé mechanické soustavě přiřadit analogický elektrický obvod a ten řešit metodami teorie obvodů, kterých je velmi mnoho a snadno se aplikují. Jen je třeba si uvědomit, že obvody jsou navzájem duální, seriovému spojení mechanických prvků přísluší paralelní zapojení analogických elekrických prvků. Detaily lze nalézt v literatuře [12]. m

Daleko větší význam má elektromechanická analogie pro analýzu buzení akustických či ultrazvukových vln a jejich šíření zejména ve vlnovodech. Vlnovod lze nahradit analogickým elekrickým obvodem s rozprostřenými parametry. Důležitý je převod mezi elektrickou a mechanickou částí systému. Tu zprostředkuje ideální elektromechanický transformátor při piezoelektrickém buzení, nebo ideální elektromechanický gyrátor pro buzení magnetostrikcí. Ten se používá také u náhradního obvodu reproduktoru. Pro piezoelektrické měniče se používá známý Masonův náhradní obvod. Podrobnosti jsou v literatuře, mj. [12]. Masonův náhradní obvod má mnoho speciálních aplikací, o jedné z nich je zmínka v monografii [1]. Do této příručky byla tato část zahrnuta zejména proto, že obor mechatronika je spojením mechaniky a elektroniky. Proto by jeho absolventi měli být s užitečnou a účinnou metodou modelování elektromechanických, ale zejména elektroelastických, soustav seznámeni. Současný výzkum se zaměřuje i do této oblasti a dosahuje nestandardní zajímavé výsledky.

6.2 Analogie mezi transportními jevy K transportním jevům dochází v případě, kdy prostorové rozložení fyzikální veličiny není homogenní a tato veličina má možnost se přenášet. Např. při nerovnoměrném rozložení teploty mají oblasti o vyšší teplotě vyšší vnitřní energii. Pak může34 docházet k přenosu tepla (tepelné energie) z míst o vyšší teplotě do míst o nižší teplotě35 . Mezi základní transportní jevy patři: • Difuze, jejíž podstatou je nerovnoměrná hustota (koncentrace), mechanismus spočívá v přenosu hmotnosti a výsledkem je tok látky. • Vnitřní tření (viskozita), které lze vyjádřit jako transport hybnosti ve směru jejího gradientu, zjednodušeně tedy gradientu rychlosti. Protože hybnost je vektorová veličina, vyjádří 33

V prvním přiblížení se předpokládá, že útlum je úměrný rychlosti. Přesně řečeno, útlum závisí i na jiných veličinách. Toto je však užitečná idealizace, podobná té, která se užívá u elektrického obvodu. 34 V případě, že není přítomen ideální tepelný izolátor. 35 Podle druhé věty termodynamické opačný přenos není možný

68

se gradient pro každou složku a tok je pak určen devíti složkami, které odpovídají tenzoru napětí. Napětí (síla) je fyzikálním významem toku hybnosti. Názorný je případ ustáleného toku v potrubí (kanále) s rychlostí v jedné ose a gradientem v příčném směru, což je doprovázeno působením tečné (třecí) síly mezi vrstvami“ kapaliny s různou rychlostí36 . ” • Vedení tepla je způsobeno přenosem energie částic v důsledku nerovnoměrného rozložení teploty. Výsledkem je tok tepelné energie, jak jsme již vysvětlili na začátku této části. • Vedení elektrického proudu, ke kterému dochází zásluhou přenosu náboje. Základní rovnice pro všechny tyto jevy má z matematichého hlediska stejný tvar, viz tab. 1. − → − → Poslední rovnice je Ohmův zákon v diferenciálním tvaru i = γ E . Znaménko minus u gradientu je nutné, poněvadž tok probíhá ve směru spádu, z místa vyšší hodnoty (počátek) do místa nižší hodnoty gradientu (konec). To jsme podrobně vysvětlili v obecné části o potenciálu 4.4.5. Vnitřní tření Difuze Vedení tepla Elektrický proud

− → q p = −Kp grad p − → q ρ = −Kρ grad ρ − → q w = −Kw grad T − → q = −γ grad φ i

Tabulka 1: Analogie popisu fyzikálních veličin. Kvůli formálnímu zjednodušení zápisu je v případě Vnitřního tření“ uvažováno pole rychlostí s nenulovou složkou pouze v 1 směru. ” Všechny tyto rovnice můžeme sjednotit do jedné obecné rovnice − → q = −K grad u.

(130)

Je zřejmé, že symbolická veličina u a jí odpovídající konkrétní veličiny mají význam potenciálu a jejich gradient pak význam intenzity. U vedení proudu se přímo jedná o intenzitu elektrického → pole. Obecný symbol − q přestavuje tok a K je koeficient. Fyzikální význam veličin v uvedených v Tabulce 1 je je popsán v Tabulce 2. Transportní jev Difuze Vnitřní tření Vedení tepla

Nehomogenita Hustota Hybnost Teplota

Přenášená veličina Hmotnost částice Hybnost částice Energie částice

Tok Hmotnost (Látka) Tečné napětí Energie

Vedení proudu

Elektrický potenciál

Náboj

Proud

Koeficient Difuze Viskozita Tepelná vodivost Vodivost

Zdroj Hmotnost Hybnost Energie Napětí

Tabulka 2: Přehled fyzikálních veličin u transportních jevů. Transportní jev Difuze Vnitřní tření Vedení tepla Vedení proudu

Nehomogenita grad ρ kg · m−4 grad p kg · s−1 grad T K · m−1 grad φ V · m−1

Přenášená veličina m kg − → p kg · m · s−1 W J = kg · m2 · s−2 Q C

qp qp qw i

Tok kg · m−2 · s−1 kg · m−1 · s−2 W · m−2 A · m−2

Kp Kp Kw γ

Koeficient m2 · s−1 kg · m−1 · s−1 W · m−1 · K −1 S · m−1

Tabulka 3: Analogie veličin, procesů a jednotek pro transportní jevy. Význam symbolů pro jednotlivé transportní jevy je v Tabulce 3. Pokud jde o koeficienty předpokládáme v prvním přiblížení, že jsou konstantní, nezávislé na žádných parametrech. Tedy jde o rovnice lineární. Ve skutečnosti jsou všechny koeficienty závislé na teplotě, případně dalších parametrech prostředí. Mohou též záviset na intenzitě budící veličiny, pak jsou vztahy nelineární. 36 Rozdíl rychlostí tedy vyvolá transport hybnosti v příčném směru mezi vrstvami kapaliny. Tím dojde k urychlení pomalejší vrstvy a zpomalení rychlejší vrstvy.

69

V závislosti na materiálu může být konstanta úměrnosti buď skalár (pro izotropní materiál), nebo tenzor (pro anizotropní materiál) - v tom případě nemusí být vektor gradientu a vektor toku rovnoběžné. Právě na příkladu elekrické vodivosti jsme v části 5.1 demonstrovali základní vlastnosti anizotropního prostředí. Průběh transportních jevů závisí též na typu prostředí. Například difuze, vedení tepla a vedení elektrického proudu existuje ve všech skupenstvích, vnitřní tření jen v plynech a kapalinách. Zejména difuze probíhá v různých skupenstvích řádově odlišnou rychlostí. V kinetické teorii plynů se na základě vlastností molekul odvozují vztahy pro příslušné koeficienty. V důsledku přenosu se rozdíly koncentrací vyrovnávají a nakonec dojde k rovnovážnému stavu s homogenní koncentrací. Relaxační doba, během níž se rozdíly vyrovnají, je různá: od řádu mikrosekund po léta. Aby se trvale udržela nehomogenita, musí v prostředí existovat vnější zdroj, který ji zajistí. Charakter zdroje závisí na typu prostředí. Pokud je prostředí nekonečné, postačují jen zřídla, případně jen nory37 . To ale není praktický případ. V praxi je prostředí vždy omezené. Pokud je prostředí omezené, musí být přítomna jak zřídla, tak nory, a musí platit podmínka rovnováhy, např. teplo dodané zřídly se musí rovnat teplu odvedenému norami, což je triviální formulace zákona zachování energie38 . Bez ohledu na charakter zdroje předpokládáme, že působí v určité oblasti a v této oblasti jej můřeme popsat symbolickou rovnicí → div − q = h,

(131)

kde h je prostorová hustota zdrojů39 . Jejich význam je v tabulce 2. Pak pro libovolný počáteční stav popsaný potenciálovým polem dojde po odeznění přechodového děje k ustálenému stavu, během něhož se již potenciálové pole s časem nemění.

6.3 Poissonova a Laplaceova rovnice Rovnicím pro transportní jevy je formálně ekvivaletní rovnice z elektrostatiky popisující vztah − → mezi elekrickým potenciálem φ a elekrickou indukcí D → − − → D = ϵ E = −ϵ grad φ,

(132)

− → kde E je intenzita elektrického pole a ϵ je absolutní permitivita40 . Formálně, podle zobecněné rovnice (130) můžeme elektrickou indukci považovat za tok, i když její jednotka C/m2 tomu neodpovídá. Z Maxwellových rovnic plyne, že zdrojem elekrické indukce je volný objemový náboj, − → div D = ρ,

(133)

kde ρ je objemová hustota volného náboje. Tato rovnice odpovídá zobecněné rovnici pro zdroje (131). Fyzikálně velmi vyznamné je spojení obou rovnic (132) a (133). Rovnici (132) dosadíme do (133). Po jednoduchých úpravách dostaneme Posissonovu rovnici ρ ∆φ = − , ϵ 37

(134)

Zde je ale problém se zákony zachování. Pro případ vedení proudu to není již na první pohled v pořádku. Pro jednoduchost uvažujme drát. Pak na jeho začátku je zřídlo náboje a na konci nora. Ve skutečnoti existuje uzavřený obvod a pohyb náboje zajišťují zdroje elektromotorického napětí, které se také ve formě prostorové hustoty zahrnují do Ohmova zákona v diferenciálním tvaru, viz např. monografii [3]. Vedení proudu má tedy poněkud jiný charakter než ostatní transportní jevy. 39 Zdrojem rozumíme zřídlo nebo noru. Bohužel to není ideální terminologie. 40 Uvažujeme tedy prostředí lineární a měkké, což je v praxi obvykle splněno. 38

70

kde ∆ je Laplaceův operátor, viz část 4.5. Pokud vyšeřujeme elektrostatické pole mimo oblast zdrojů, přejde Poissonova rovnice v rovnici Laplaceovu ∆φ = 0.

(135)

Poissonovu a Laplaceovu rovnici lze úplně stejně odvodit pro všechny transportní jevy. Transportní jevy v ustáleném stavu lze tedy řešit úplně stejnými metodami jako jevy elektrostatické. Tato analogie má obrovský význam v tom, že v oblasti elektrostatiky byla intenzivním výzkumem v minulosti vypracována řada sofistikovaných metod, jak řešit Poissonovu a Laplaceovu rovnici. Ty lze okamžitě převzít pro efektivní řešení transportních jevů. Existují však i chytře vymyšlené přímé metody jak elektrostatické pole řešit ve speciálních případech. Pokud geometrie pro transportní jev odpovídá geometrii v elektrostatice, lze tuto metodu rovněž převzít. Tato možnost zůstává většinou nevyužita. Pro vektorový potenciál lze rovněž odvodit Poissonovu a Laplaceovu rovnici. Pro případ magnetického pole lze ukázat [3], že má tvar → − − → ∆ A = −µ i ,

(136)

− → kde i je objemová hustota volných proudů a µ je absolutní permabilita41 . Pro případ izotropního − → prostředí se jedná o tři rovnice pro složky vektoru A, každá z nich je analogická rovnici pro skalární potenciál. I zde lze tedy v širokém rozsahu aplikovat metody elektrostatiky.

41

Opět uvažujeme prostředí lineární a měkké, i když to pro feromagnetické látky neplatí.

71

72

7 Postupy řešení matematických modelů V předcházejících kapitolách byla uvedena řada druhů matematického popisu fyzikálních jevů. Tato kapitola představuje vybrané metody a nástroje, kterými je možné správně zformulované matematické modely řešit a získat tím z modelů požadované informace. Nejprve je v kapitole krátce shrnuto, jaké typy popisu je možné a vhodné použít pro formulaci modelů různých druhů přírodních, technologických nebo výrobních procesů. Dále jsou zmíněny druhy soustav rovnic, které je často nutno využít pro popis vývoje složitějších systémů.

7.1 Klasifikace rovnic a jejich použití Na začátku si vytyčme předpoklad, že většina přírodních, fyzikálních, chemických a dalších procesů směřuje k ustavení určité, statické nebo dynamické rovnováhy v systému. Různé druhy dějů je zvykem, v závislosti na rychlosti jejich průběhu, popisovat různými druhy matematických rovnic. Rychlé děje, u kterých není nutné a nebo je příliš složité sledovat jejich průběh, jsou často popisovány algebraickými rovnicemi42 . Děje, u kterých je žádoucí pozorovat, popisovat a předpovídat jejich průběh, bývají popisovány buď obyčejnými diferenciálními nebo parciálními diferenciálními rovnicemi. Popis stejného jevu jedním z uvedených typů rovnic se může případ od případu lišit podle toho jaká informace je předmětem modelářova zájmu. Velkou část přírodních dějů lze popsat poměrně jednoduchými lineárními parciálními diferenciálními rovnicemi. V obecném případě v nich vystupují čtyři nezávisle proměnné, tři souřadnice a čas. Analytické řešení parciálních diferenciálních rovnnic je, i při redukci počtu nezávisle proměnných, možné jen pro velmi jednoduché případy, které nicméně dávají klíčové fyzikální výsledky, nebo naleznou praktické uplatnění. K dosažení větší názornosti jsou v dalším textu často voleny co nejjednodušší modelové příklady z teorie obvodů. Ta byla zformulována pro velmi důležitou skupinu kvazistacionárních elektromagnetických jevů. Tyto jevy jsou řešeny s velkou přesností a mnohem jednoduššími prostředky než by vyžadovala obecná nestacionární teorie pole. Při vyšších frekvencích zdroje lze využít modely obvodů s rozloženými parametry, které jsou při harmonickém buzení popsány parciálními diferenciálními rovnicemi. Ve zjednodušování modelů můžeme jít ještě dále. Často dokážeme všechny podstatné informace o poměrně složitých obvodech získat řešením soustavy lineárních algebraických rovnic. 7.1.1 Algebraické a transcendentní rovnice Skoro celá středoškolská fyzika používá k popisu přírodních jevů buď algebraické (polynomiální) nebo transcendentní (nepolynomiální) rovnice. Jde vlastně o jednoduché matematické modely těchto jevů. Algebraické rovnice možňují popsat aktuální stav simulovaného systému, charakteristiky systému, popřípadě ustálený stav. Ilustrace využití transcendentních a algebraických rovnic bude dále demonstrována na příkladu popisu obvodu střídavého proudu prostřednicvím Komplexní symbolické metody“. ” Komplexní symbolická metoda Mějme elektrický obvod buzený harmonickým napětím nebo proudem. Harmonické buzení můžeme popsat několika způsoby, např.: u(t) = Uo sin ωt,

u(t) = Uo cos ωt,

i(t) = Io sin ωt,

i(t) = Io cos ωt,

(137)

42 Volba typu rovnice podle rychlosti změn ve sledovaném systému není striktním pravidlem. Například kvazistacionární děje jsou děje relativně pomalé. Jde obvody se soustředěnými parametry, prakticky do 10 MHz. Přitom k jejich popisu stačí algebraické rovnice. Naopak pro rychlejší děje (obvody s rozloženými parametry) se používají parciální diferenciální rovnice pro dvě proměnné, souřadnici a čas. Pro harmonické buzení pak postačuje obyčejná diferenciální rovnice pro souřadnici.

73

kde ω = 2πf je úhlová (kruhová) frekvence, Uo je (maximální) amplituda napětí a Io je (maximální) amplituda proudu. Pro jednoduchost uvažujeme fázovou konstantu rovnou nule. K pozdějšímu použití si nyní zavedeme popis pasivních dvojpólů prostřednictvím funkčních předpisů mezi budícím proudem i(t) jako podnětem a svorkovým napětím u(t) jako odezvou na podnět. 1. Pro resistor s odporem R platí u(t) = Ri(t).

(138)

2. Pro induktor s indukčností L platí u(t) = L

di(t) . dt

(139)

3. Pro kapacitor s kapacitou C platí 1 u(t) = C

∫

t

i(τ )dτ.

(140)

t=0

Nevýhodou zápisu (137) je skutečnost, že po dosazení např. třetího výrazu z (137) do výše uvedeného vztahu (139) pro napěťovou odezvu na induktoru, dostaneme časovou závislost napětí jako funkci cos. Abychom mohli buzení a odezvu vzájemně porovnat, musíme tento kosinus převést na sinus přidáním fázového posuvu −π/2. V uvažovaném případě je to jednoduché, ale už u mírně složitějších obvodů může být porovnání matematicky dosti pracné. Proto je výhodné použít budící funkci jejíž derivace se liší jen koeficientem. Tuto vlastnost má exponenciální funkce. V teorii funkcí komplexní proměnné se využívá definice ejx = e(jx) = cos x + j sin x, (141) kde j je komplexní jednotka43 . Pro obecný harmonický průběh napětí proto platí vztah ˆo ejωt = Uo cos(ωt + φu ) + j Uo sin(ωt + φu ). Uo ej(ωt+φu ) = U kde ˆo = Uo ejφu U

(142)

je komplexní amplituda neboli fázor. Veličina u(t) = Uo ej(ωt+φu ) = Uo cos(ωt + φu ) + jUo sin(ωt + φu ).

(143)

se nazývá rotující fázor. V komplexní amplitudě jsou podle vztahu (142) obsaženy všechny informace o harmonickém průběhu: amplituda Uo i fázový posuv φu . Z rotujícího fázoru (143) získáme časový průběh napětí jako jeho reálnou nebo imaginární část. Stejným způsobem jako v případě napětí se dá zavést rotující fázor proudu i(t) = Iˆ0 ejωt .

(144)

Zde je důležité zdůraznit, že fázor nemá fyzikální význam, je to komplexní veličina. Jak uvidíme dále, je to velmi důležitá pomocná matematická veličina. Fyzikální význam mají parametry fázoru (viz (142)) stejně jako jeho reálná či imaginární část (viz (143)). 43 Nepoužíváme v matematice zavedený symbol i, protože tato značka i se obvykle používá pro okamžitou hodnotu proudu.

74

Dosadíme-li rotující fázory (143) a (144) do definičních vzorců předepisujících napěťové odezvy na budící proudy (138) až (140), po jednoduchých úpravách získáme následující vztahy: 1. pro resistor ˆo = RIˆo = ZˆR Iˆo , U

kde

ZˆR = R,

(145)

ˆo = jωLIˆo = ZˆL Iˆo , U

kde

ZˆL = jωL,

(146)

1 ˆ Io = ZˆC Iˆo , jωC

kde

ZˆC =

2. pro induktor 3. pro kapacitor ˆo = U

1 . jωC

(147)

Veličiny ZˆR , ZˆL , ZˆC se nazývají po řadě impedance rezistoru, induktoru a kapacitoru. Význam symbolických rovnic (145) až (147) vysvitne po dosazení fázoru proudu Iˆo . Předpokládáme, že jeho fázová konstanta je nulová, pak Iˆo = Io . Pro fázor napětí na rezistoru dostaneme podle (145) vztah ˆo = ZˆR Iˆo = RIo . U Fázor napětí je reálný, fázový posuv je nulový, napětí je ve fázi s proudem. Amplituda napětí je úměrná odporu. Pro fázor napětí na induktoru dostaneme ze (146) vztah44 ˆo = ZˆL Iˆo = jωLIo = ωLejπ/2 Io . U Z toho plyne, že napětí na induktoru předbíhá proud o π/2 = 90o . Amplituda napětí je úměrná frekvenci a indukčnosti. Pro fázor napětí na kapacitoru dostaneme ze (147) vztah45 ˆo = ZˆC Iˆo = U

1 1 −jπ/2 Io = e Io . jωC ωC

Odtud vidíme, že napětí na kapacitoru se zpožďuje za proudem o úhel π/2 = 90o . Amplituda napětí je nepřímo úměrná frekvenci a kapacitě. Vztahy (145) až (147) jsou speciálními případy Ohmova zákona pro střídavé proudy ˆo U Iˆo = , Zˆ

(148)

kde Zˆ je impedance. Získá se z impedancí elementárních prvků stejně jako v případě rezistorů. Důležité je si uvědomit, že fázor by měl být v exponenciálním tvaru. Pak se z exponentu určí fázový posuv v radiánech. Amplituda fázoru je číslo před exponentem. Ve složkovém tvaru je to absolutní hodnota. Fázový posuv pak lze určit z reálné a imaginární složky. Komplexní symbolická metoda platí jen pro lineární kombinace. Nelze ji přímo použít na výpočet výkonu, poněvadž zde je součin exponenciál. Naštěstí lze použít vztah ˆo Iˆo∗ Pˆ = U

(149)

kde Iˆo∗ je komplexně sdružený fázor proudu a Pˆ je komplexní výkon. Jeho reálná část je činný výkon a imaginární část je jalový výkon. Zdánlivý výkon a účiník pak vypočteme z těchto výkonů. Pro jednoduché obvody se dá symbolický počet aplikovat přímo. Aplikace symbolického počtu vede k získání analytických výsledků. Níže je uveden příklad obvodu vhodného pro analytické řešení. 44 45

Použijeme vztahu j = ejπ/2 , který plyne např. ze vztahu (141). Použijeme vztahu 1j = −j = e−jπ/2 , který plyne např. ze vztahu (141).

75

PŘÍKLAD 23 Mějme seriový RLC obvod z Kapitoly 6 zobrazený pro připomenutí na obrázku 32. Uvažovaný obvod je buzený harmonickým napětím o amplitudě Uo , frekvenci ω a nulové fázové konstantě. Modelový případ obvodu popíšeme pomocí Ohmova zákona pro střídavé proudy (148). Impedance jednotlivých prvků je dána vztahy (145) až (147). Výsledná impedance je jejich součtem. ( ) 1 1 ˆ ˆ ˆ ˆ Z = ZR + ZL + ZC = R + jωL + = R + j ωL − (150) jωC ωC

Obrázek 32: Sériový RLC obvod. Absolutní hodnota (velikost) impedance je √ ) ( 1 2 2 Z = R + ωL − ωC

(151)

Pro fázový úhel výsledné impedance dostaneme ze vztahu (150) výraz tan φ =

ωL − R

1 ωC

.

Z posledního výrazu automaticky vyplývá, že proud tekoucí obvodem může být vůči budícímu napětí posunut v rozmezí od −90◦ do +90◦ v závislosti na frekvenci a hodnotách prvků. Pro proud tekoucí obvodem dostáváme jednoduchý výraz U ˆ o Uo Io = |Iˆo | = = Zˆ Z Jakmile známe fázor proudu, můžeme vypočítat fázory napětí (tj. amplitudu a fázový posuv) na jednotlivých prvcích obdobným způsobem. Pomocí vztahu (149) lze pak vypočítat odebíraný výkon a jeho parametry, případně výkony na jednotlivých prvcích. Ty ale budou buď jen činné nebo jen jalové. Pro zvolený modelový příklad se jeví demonstrovaný popis jako dostatečný. V případech složitějších obvodů a přírodních procesů však bude třeba využít jiné matematické nástroje, které umožní použití výpočetní techniky namísto lidské pracovní síly, čímž dojde k mnohonásobnému zrychlení práce s modely.

76

7.1.2 Obyčejné diferenciální rovnice Jestliže je předmětem modelování vývoj vybrané veličiny v čase, pak se v matematickém popisu jevu vyskytnou s velkou pravděpodobností obyčejné diferenciální rovnice. V obyčejných diferenciálních rovnicích figurují členy s derivacemi hledané veličiny podle jedné nezávisle proměnné (často podle času). Diferenciální rovnice mohou být různého řádu. Řád rovnice je dán nejvyšším řádem derivace, která v rovnici vystupuje. Níže jsou uvedeny příklady obyčejných diferenciálních rovnic prvního a druhého řádu.

PŘÍKLAD 24 Uvažujme dva jednoduché obvody střídavého proudu. V prvním z nich je zapojen rezistor v sérii s cívkou a ve druhém potom rezistor v sérii s kondenzátorem. Oba dva obvody jsou níže popsány obyčejnými diferenciálními rovnicemi prvního řádu.

(a) RL obvod.

(b) RC obvod.

Obrázek 33: Obvody popsané obyčejnými diferenciálními rovnicemi prvního řádu.

V případě RC obvodu (viz obrázek 33(b)) popsaného obyčejnou diferenciální rovnicí prvního řádu je hledanou veličinou elektrický náboj q(t).

V případě RL obvodu (viz obrázek 33(a)) popsaného obyčejnou diferenciální rovnicí prvního řádu je hledanou veličinou elektrický proud i(t). L

di + Ri = U0 sin(ωt + φ0 ) dt

R

1 dq + q = U0 sin(ωt + φ0 ) dt C

PŘÍKLAD 25 Využití obyčejných diferenciálních rovnic druhého řádu je níže ilustrováno na příkladech RL a RLC obvodů střídavého proudu (viz obrázek 34). V obou případech je tentokrát hledanou veličinou elektrický náboj q(t).

77

(a) RL obvod.

(b) RLC obvod.

Obrázek 34: Obvody popsané obyčejnými diferenciálními rovnicemi druhého řádu.

L

d2 q dq = U0 sin(ωt + φ0 ) +R 2 dt dt

L

d2 q 1 dq + R + q = U0 sin(ωt + φ0 ) 2 dt dt C

Ještě jednodušší příklady obyčejných diferenciálních rovnic lze nalézt v mechanice, kde rychlost je derivací polohy podle času a zrychlení je časovou derivací rychlosti. 7.1.3 Parciální diferenciální rovnice Při konstrukci složitějších modelů obvykle nestačí popisovat časový vývoj veličiny na jednom izolovaném“ místě. Jinými slovy, do modelu je třeba zavést další nezávisle proměnné, prostorové ” souřadnice, jejichž pomocí je zohledněno rozložení zkoumané veličiny v blízkém okolí sledovaného místa. Toto rozšíření modelu vede v matematické formulaci na popis prostřednictvím parciálních diferenciálních rovnic. V modelování často používané parciální diferenciální rovnice druhého řádu dělíme do tří typů, z nichž všechny mají konkrétní fyzikální aplikace. eliptická parabolická hyperbolická

∂2u ∂2u + 2 = f, ∂x2 ∂y ∂u ∂ 2 u K + 2 = f, ∂t ∂x ∂2u ∂2u K 2 − 2 = f. ∂t ∂x

(152) (153) (154)

Výše jsou zapsány jen speciální jednoduché tvary, ale zařadit do jednoho z uvedených typů lze každou lineární parciální diferenciální rovnici 2. řádu (např. se smíšenými druhými a prvními derivacemi). Symbolem u je označena hledaná proměnná (vyšetřovaná veličina), často to bývá potenciál. Konstanta K je odvozena z příslušných výchozích rovnic46 . Funkce f na pravé straně má různý význam podle druhu řešeného problému. U transportních jevů to může být objemová hustota zdrojů, v elektrostatice hustota objemového náboje, v magnetizmu objemová hustota proudu. Pro jednoduchost je uvedena parciální derivace jen podle jedné či dvou souřadnic. Obvykle je parciální derivace provedena podle všech tří souřadnic. Derivace podle souřadnic jsou druhého 46 U hyperbolické rovnice se z ní může odvodit rychlost šíření vlny. U parabolické rovnice má význam převrácené hodnoty koeficientu vedení tepla nebo koeficintu difúze.

78

řádu. Rovnice se liší derivací podle času. Eliptické rovnice časovou derivaci neobsahují. Parabolické rovnice mají první derivaci podle času a hyperbolické rovnice obsahují druhou derivaci podle času. Zařazení rovnic v rámci této klasifikace má význam pro správnou formulaci okrajových (a počátečních) podmínek a umožňuje dále využít analogie a vzájemné souvislosti. Z matematického hlediska má každý ze tří uvedených typů parciálních diferenciálních rovnic nekonečně mnoho navzájem odlišných analytických řešení. K jednoznačnému řešení je nutno dodat okrajové a v případě parabolických či hyperbolických rovnic i počáteční podmínky (viz kapitola 7.3.4). Eliptické parciální diferenciální rovnice Eliptické parciální diferenciální rovnice popisují stacionární (tj. na čase nezávislé) fyzikální jevy.

PŘÍKLAD 26 Příkladem eliptické parciální diferenciální rovnice může být Poissonova nebo Laplaceova rovnice. Poissonova rovnice představuje popis elektrického pole ve 2D oblasti se zdroji (volný elektrický náboj). ∆φ =

Laplaceova rovnice popisuje elektrického pole ve 2D oblasti bez zdrojů. ∆φ =

∂2φ ∂2φ ρ + =− 2 2 ∂x ∂y ϵ

∂2φ ∂2φ + =0 ∂x2 ∂y 2

Na eliptickou parciální diferenciální rovnici přecházejí ostatní rovnice tehdy, když simulovaná veličina dosáhne v modelu ustáleného stavu. Parabolické parciální diferenciální rovnice Parabolické rovnice popisují nevratné, nestacionární jevy. Parciální diferenciální rovnice s časovými členy vycházejí zpravidla z názorných integrálních forem zákonů zachování hmotnosti, energie, náboje, atp.

PŘÍKLAD 27 Typickými praktickými příklady jevů popisovaných parabolickými parciálními diferenciálními rovnicemi jsou difúze a vedení tepla. Pro tyto případy nevykazují řešení žádnou periodicitu a jsou stabilní. Analogií k matematickému popisu difúze je model vedení tepla. Vedení tepla je možné pozorovat prostřednictvím sledování změn teploty T . K toku − → q w tepla dochází proti směru gradientu teploty ∇T . Parametr KT reprezentuje materiálové vlastnosti prostředí.

Parabolickou diferenciální rovnici lze použít k popisu vyrovnávání hustoty roztoku ρ prostřednictvím difúze. Časová změna hustoty ∂ρ ∂t roztoku ve sledovaném objemu odpovídá celkové prostorové změně difuzního toku hus→ toty − q ρ tamtéž. Parametr Kρ reprezentuje materiálové vlastnosti prostředí.

79

− → q ρ = −Kρ ∇ρ ∂ρ ∂t + Kρ ∇ ( · ∇ρ = 0 ∂ρ ∂t

+ Kρ

∂2ρ ∂x2

+

∂2ρ ∂y 2

+

∂2ρ ∂z 2

− → q w = −KT ∇T ∂T ∂t + KT ∇ ( · ∇T = 0

)

∂T ∂t

=0

+ KT

∂2T ∂x2

+

∂2T ∂y 2

+

∂2T ∂z 2

) =0

Hyperbolické parciální diferenciální rovnice Hyperbolické rovnice popisují vlnění. V elektromagnetizmu to jsou známé netlumené nebo jen slabě tlumené elektromagnetické vlny. V elasticitě to jsou elastické (akustické, ultrazvukové) vlny.

PŘÍKLAD 28 Pro netlumené elektromagnetické vlnění ve vakuu a za nepřítomnosti zdrojů lze z Maxwellových rovnic odvodit vztahy ∂2E ∂t2 ∂2B ∂t2

Šíření zvukové vlny prostorem je možné popsat popsat rovnicí ∆p −

− c2 ∆E = 0 − c2 ∆B = 0

1 ∂2p =0 c2 ∂t2

ve které p reprezentuje akustický tlak a malé c představuje rychlost šíření vlny prostorem.

kde c značí rychlost světla ve vakuu.

7.1.4 Diferenciálně algebraické rovnice/soustavy Speciální třída rovnic, resp. soustav rovnic nese označení diferenciálně algebraické rovnice (DAE - Differential Algebraic Equation). Jde vlastně o diferenciální rovnice doplněné určitou speciální podmínkou, která připouští jen vybraná z množiny všech řešení. Podmínka/y vymezující množinu přípustných řešení mají v případě DAE podobu buď algebraického výrazu nebo algebraické rovnice. Diferenciálně algebraické rovnice jsou jakýmsi dalším krokem zobecňování diferenciálních rovnic ve směru: homogenní → nehomogenní → diferenciálně-algebraické. Popisovaný typ rovnic/soustav se používá například při modelování elektrických obvodů s nelineárními prvky. Častým příkladem využití DAE v mechanice je popis pohybu kyvadla u kterého je trajektorie pohybu vymezena délkou jeho ramene. Konkrétní příklady ani způsoby řešení diferenciálně-algebraických rovnic nebudou v této knize uváděny. Případní zájemci mohou nalézt příklady a implementaci modelů vybraných jevů pomocí softwarového balíku Sundials (viz [4]) na adrese https://computation.llnl.gov/casc/sundials/documentation/stb_guide/ sundialsTB.html.

7.2 Analytické řešení rovnic a soustav rovnic V případě, že potřebujeme matematicky popsat systém, jehož vývoj/stav je ovlivněn vícero vzájemně nezávislými procesy, pak vede tento popis často na soustavu rovnic jednoho z následujících typů. • AE (Algebraic Equations) - algebraické rovnice • ODE (Ordinary Differential Equations) - obyčejné diferenciální rovnice 80

• PDE (Partial Differential Equations) - parciální diferenciální rovnice • DAE (Differential Algebraic Equations) - diferenciálně-algebraické rovnice Soustavu popisující uvažovaný systém můžeme řešit různými způsoby, které volíme na základě druhu použitých rovnic, časové náročnosti řešení soustavy a podle dostupných programových nástrojů. Tradičním postupem získání hledané informace je analytické řešení problému. Jeho výsledkem je často funkce, vzoreček, popisující vývoj vybrané veličiny v čase a/nebo v prostoru. Analytické řešení je výhodné z hlediska své názornosti a srozumitelnosti, ale lze ho nalézt jen pro jednoduché případy modelů. Níže jsou představeny způsoby řešení soustav algebraických rovnic, soustav obyčejných diferenciálních rovnic a řešení parciálních diferenciálních rovnic. V rámci této kapitoly jsou uvedeny ilustrativní příklady způsobů analytického řešení jednoduchých matematických modelů. Příklady jsou řešeny jednak postupy, které bychom použili při počítání na papíře a jednak za pomoci simulačního nástroje Sage (viz [19]). Alternativou k Sage může být komerční software Maple (viz [13]). 7.2.1 Analytické řešení soustav algebraických rovnic Modely vzájemně provázaných dějů vedou často buď přímo nebo po určitých úpravách k soustavě algebraických (polynomiálních) rovnic. Velmi významná část modelů je tvořena soustavami lineárních algebraických rovnic. Numerické modely řešené na počítači vyžadují obvykle formulaci úlohy právě ve tvaru soustavy lineárních algebraických rovnic. Formulace modelu ve tvaru soustavy lineárních algebraických rovnic může být získána buď na základě lineární povahy rovnic popisujících simulovaný děj nebo pomocí linearizace původně nelineárního matematického popisu. Převod parciální diferenciální rovnice na soustavu lineárních algebraických rovnic lze provést například prostřednictvím diskretizace úlohy v čase a v prostoru. Následuje příklad soustavy lineárních algebraických rovnic popisujicích proudy a napětí v jednoduchém obvodu se soustředěnými parametry.

PŘÍKLAD 29 Mějme jednoduchý elektrický obvod sestavený ze dvou zdrojů a ze tří rezistorů (viz obrázek 35).

Obrázek 35: Modelový obvod. Soustava pro výpočet proudů I1 , I2 , I3 tekoucích obvodem má tvar −I1 + I2 − I3 = 0, R1 I1 + R2 I2 = U1 , R2 I2 + R3 I3 = U2 . Maticový zápis soustavy rovnic, která popisuje velikosti proudů v obvodu je   −1 1 −1 0  R 1 R 2 0 U1  . 0 R 2 R 3 U2 81

Standardním postupem analytického řešení soustav lineárních algebraických rovnic je Gaussova eliminace. Řešení soustavy s využitím Gaussovy eliminace zde nebude provedeno. Namísto toho je níže uvedený krátký úsek kódu pro nalezení řešení v programu Sage. sage: (I_1, sage: sage: sage: sage: sage:

var('I_1 I_2 I_3 R_1 R_2 R_3 U_1 U_2') I_2, I_3, R_1, R_2, R_3, U_1, U_2) eq1 = (-1)*I_1 + I_2 - I_3 == 0 eq2 = R_1*I_1 + R_2*I_2 == U_1 eq3 = R_2*I_2 + R_3*I_3 == U_2 proudy = solve([eq1, eq2, eq3], I_1,I_2,I_3) print proudy

Analytické řešení uvažované soustavy je ve výše uvedeném kódu uloženo do proměnné proudy“ ” a vychází jako

I1 =

(R2 + R3 )U1 − R2 U2 R3 U1 + R1 U2 R 2 U1 − R 1 U2 − R 2 U2 , I2 = , I3 = − . (R1 + R3 )R2 + R1 R3 (R1 + R3 )R2 + R1 R3 (R1 + R3 )R2 + R1 R3

V případě se zadanými hodnotami napětí zdrojů a velikostí odporů je třeba definovat předepsané hodnoty jako parametry (přiřazením hodnoty) a ne jako proměnné (var(′ . . .′ )). sage: var('I_1 I_2 I_3') (I_1, I_2, I_3) sage: R_1 = 50.0 sage: R_2 = 75.0 sage: R_3 = 100.0 sage: U_1 = 12.0 sage: U_2 = 24.0 sage: eq1 = (-1)*I_1 + I_2 - I_3 == 0.0 sage: eq2 = R_1*I_1 + R_2*I_2 == U_1 sage: eq3 = R_2*I_2 + R_3*I_3 == U_2 sage: proudy = solve([eq1, eq2, eq3], I_1,I_2,I_3) sage: print proudy [ [I_1 == (6/325), I_2 == (48/325), I_3 == (42/325)] ]

Proudy tekoucí jednotlivými větvemi modelového obvodu mají velikosti [I1 ; I2 ; I3 ] = [0, 0184615384615385; 0, 147692307692308; 0, 129230769230769] A. Na soustavě z uvedeného příkladu budou v dalším textu předvedeny vybrané metody numerického řešení soustav lineárních algebraických rovnic.

7.2.2 Analytické řešení soustav obyčejných diferenciálních rovnic Na soustavu obyčejných diferenciálních rovnic vede, mimo jiné, řešení diferenciálních rovnic vyšších řádů.

PŘÍKLAD 30 Příklady obyčejných diferenciálních rovnic vyššího, konkrétně druhého, řádu jsou rovnice kmitání mechanického (viz obr. 36(b)) nebo elektromagnetického (viz obr. 36(a)) oscilátoru. Níže jsou obě analogické rovnice uvedeny. Je uvažován případ vlastního, ne nuceného, kmitání. 82

(a) LC oscilátor.

(b) Mechanický oscilátor.

Obrázek 36: Elektromechanická analogie. d2 x + kx = 0 dt2 Obdobně jako u případu elektromagnetického oscilátoru vyřešíme příslušnou charakteristickou rovnici

d2 q 1 + q=0 2 dt C Rovnice s odpovídajícím charakteristickým polynomem a její řešení následuje. 1 = λ2 + LC √0

m

L

λ1,2 = ±j

λ2 +

1 LC

λ1,2

Z kořenů chararakteristické rovnice plyne řešení diferenciální rovnice popisující kmitání LC oscilátoru. Protože diferenciální rovnici musí splňovat reálná i imaginární složka řešení

k m

=√ 0 k = ±j m

a stanovíme řešení obyčejné diferenciální rovnice. (√

x(t) = k1 cos

ej(λt) = cos(λt) + j sin(λt),

k t m

)

(√

+ k2 sin

k t m

)

Pro jednoznačnost řešení je třeba zohlednit počáteční výchylku x(0) a její derivaci x′ (0).

má toto (obecné) řešení tvar t t q(t) = k1 cos( √ ) + k2 sin( √ ), LC LC kde {k1 , k2 } jsou reálné parametry. Aby bylo řešení jednoznačné, je třeba doplnit ještě počáteční podmínky q(0), q ′ (0).

Demonstrace řešení diferenciální rovnice popisující oscilace v LC obvodu je níže provedena s využitím programu Sage. sage: var('L C t') (L, C, t) sage: q = function('q',t) sage: de = L*diff(q,t,2) + q/C == 0 sage: desolve(de, q, ivar=t) --------------------------------------------------------------------------TypeError Traceback (most recent call last)

83

. . . TypeError: ECL says: Maxima asks: Is C*L positive or negative? --------------------------------------------------------------------------sage: (C > 0).assume() sage: (L > 0).assume() sage: desolve(de, q, ivar=t) k2*cos(t/(sqrt(C)*sqrt(L))) + k1*sin(t/(sqrt(C)*sqrt(L)))

Výsledky získané za použití lidké síly se shodují s výsledky z programového řešiče. Chybová hláška je ve výpisu programu ponechána záměrně aby bylo jasné, že pokud si program řekne“ ” o doplňující informace, tak mu je musíme poskytnout.

PŘÍKLAD 31 Ke znatelnému zvýšení obtížnosti postupu řešení dojde, když do obvodu z předcházejícího příkladu přidáme zdroj střídavého napětí. Z kmitání na vlastní frekvenci se tím stane kmitání nucené. Schema upraveného obvodu ukazuje obrázek 37.

Obrázek 37: Nucené kmitání LC oscilátoru. V diferenciální rovnici se přidání zdroje projeví objevením nenulového členu na pravé straně rovnice. Tím vznikne takzvaná nehomogenní diferenciální rovnice. d2 q 1 + q = U0 sin(ωt + φ0 ) 2 dt C Nehomogenní diferenciální rovnice je řešena v několika krocích. Nejprve je řešena charakteristická rovnice příslušející diferenciální rovnici s nulovou pravou stranou. L

Lλ2 + λ1,2 Řešením homogenní rovnice je funkce (√ qH (t) = c1 sin

1 C

= √0 1 = ±j LC

) (√ ) 1 1 t + c2 cos t , LC LC

kde parametry c1 , c2 jsou reálná čísla. V dalším kroku je třeba nalézt část řešení zohledňující nenulovost pravé strany (partikulární řešení), tj. připojení zdroje. Kvůli tomu je provedena takzvaná variace konstant (variation of 84

parameters). Postup variace konstant je založen na nahrazení parametrů c1 , c2 v předpisu qH (t) funkcemi c1 (t), c2 (t). Funkce qH (t) je nejprve jednou zderivována. √ √ (√ ) (√ ) 1 1 1 1 ′ qH (t) = c1 (t) LC cos t − c (t) sin t + ... 2 LC LC (√ ) (√ )LC 1 1 ′ +c′1 (t) sin LC t + c2 (t) cos LC t

Následně je stanovena/zvolena dodatečná podmínka kladená na ci (t), i ∈ {1, 2}. Zavedení dodatečné podmínky, přidání druhé rovnice, je potřeba kvůli jednoznačnému určení hledaných funkcí ci (t), i ∈ {1, 2}. Pro výpočet dvou neznámých funkcí potřebujeme dvě rovnice. Tvar níže vybrané podmínky (155) zjednoduší strukturu výsledné soustavy a usnadní nalezení řešení. (√ 0=

c′1 (t) sin

1 t LC

(√

) +

c′2 (t) cos

1 t LC

) (155)

Po zavedení dodatečné podmínky je funkce qH (t) zderivována podruhé. (√ ) (√ ) 1 1 1 1 ′′ qH (t) = −c1 (t) LC sin t − c (t) cos t + ... 2 LC √ √ LC (√ LC) (√ ) 1 1 1 1 ′ +c′1 (t) LC cos LC t − c2 (t) LC sin LC t ′ (t), q ′′ (t) dosazeny do řešené diV následujícím kroku jsou funkce qH (t) a její derivace qH H ferenciální rovnice, která pak spolu s podmínkou (155) vytvoří soustavu rovnic pro nalezení partikulárního řešení. √ √ (√ ) (√ ) 1 1 1 1 ′ (t)L c′1 (t)L LC cos t − c sin t = U0 sin (ωt + φ0 ) 2 LC LC ) LC (√ ) (√ 1 1 ′ c′1 (t) sin LC t + c2 (t) cos LC t = 0

Následuje nalezení funkcí c1 (t), c2 (t). c′1 (t) √

=

−c′2 (t)

) (√ 1 t cos LC (√ ) 1 sin t LC

−c′2 (t)

(√ ) 1 cos2 t LC (√ ) 1 sin t LC

c′2 (t)

(√ ) 1 sin2 t LC (√ ) 1 sin t LC

− U0 sin (ωt + φ0 ) = √ (√ )] ∫ ′ ∫[ 1 c2 (t)dt = − C U sin (ωt + φ ) sin t dt = 0 0 L √ (√ ) (√ LC ) ′ 1 1 u = LC u = − sin LC t 1 cos LC t = = v = sin (ωt + φ0 ) v ′ = ω cos (ωt + φ0 ) ) [ (√ ) ] (√ ∫ 1 1 t sin (ωt + φ ) − ωU C cos t cos (ωt + φ ) dt = = U0 C cos 0 0 0 LC LC √ (√ ) (√ ) ′ 1 1 u = LC u = cos LC t 1 sin LC t = = v = cos (ωt + φ0 ) v ′ = −ω sin (ωt + φ0 ) √ (√ ) (√ ) 1 LC 1 = U0 C cos t sin (ωt + φ ) − ωU C sin t cos (ωt + φ0 ) − 0 0 LC LC √ [ (√ ) ] 1 ∫ 1 sin −ω 2 U0 C LC 1 LC t sin (ωt + φ0 ) dt C L

Argument posledního z uvedených integrálů je stejný jako argument integrálu který počítáme. To znamená, že v této chvíli můžeme vyjádřit funkci c2 (t). (√ ) (√ ) √ U0 C U0 Cω LC 1 1 c2 (t) = cos t sin (ωt + φ0 ) − sin t cos (ωt + φ0 ) 1 − ω 2 LC LC 1 − ω 2 LC LC 85

Postupem analogickým k výpočtu c2 (t) nalezneme c1 (t). c′2 (t)

=

−c′1 (t)

(√ ) 1 sin t LC (√ ) 1 cos t LC

√ (√ )] ∫[ 1 c′1 (t)dt = C U sin (ωt + φ ) cos t dt = 0 0 L √ (√ ) (√ LC ) ′ 1 1 u = LC u = cos LC t 1 sin LC t = = v = sin (ωt + φ0 ) v ′ = ω cos (ωt + φ0 ) √ (√ ) [ ( ) ] ∫ 1 1 = U0 C sin t sin (ωt + φ ) − ωU C sin t cos (ωt + φ ) dt = 0 0 0 LC √ (√ ) (√ ) LC ′ 1 1 u = LC u = − sin LC t 1 cos LC t = = v = cos (ωt + φ0 ) v ′ = −ω sin (ωt + √φ0 ) (√ ) (√ ) 1 LC 1 = U0 C sin cos LC t sin (ωt + φ0 ) + ωU0 C LC t cos (ωt + φ0 ) + √ [ (√ ) ]1 ∫ 1 cos +ω 2 U0 C LC 1 LC t sin (ωt + φ0 ) dt ∫

Z posledního uvedeného vztahu už je možné určit funkci c1 (t). (√ ) (√ ) √ U0 C U0 Cω LC 1 1 c1 (t) = sin t sin (ωt + φ0 ) + cos t cos (ωt + φ0 ) 1 − ω 2 LC LC 1 − ω 2 LC LC Partikulární řešení má tvar (√ qP (t) = c1 (t) sin

1 t LC

)

(√ + c2 (t) cos

1 t LC

) =

U0 C sin (ωt + φ0 ) 1 − ω 2 LC

Řešení obyčejné diferenciální rovnice druhého řádu popisující nucené kmitání LC-oscilátoru je součtem vypočteného obecného a partikulárního řešení. (√ ) (√ ) 1 1 U0 C q(t) = qH (t) + qP (t) = c1 sin t + c2 cos t + sin (ωt + φ0 ) LC LC 1 − ω 2 LC Výše odvozené řešení nehomogenní obyčejné diferenciální rovnice lze s vynaložením nesrovnatelně menšího úsilí získat v programu Sage příkazy zapsanými níže47 . sage: var('L C omg fi t U_0') (L, C, omg, fi, t, U_0) sage: q = function('q',t) sage: deq = L*diff(q, t, 2) + q/C == U_0*sin(omg*t + fi) sage: assume(C>0) sage: assume(L>0) sage: assume(U_0>0) sage: desolve(deq, q, ivar=t) k2*cos(t/(sqrt(C)*sqrt(L))) - C*U_0*sin(omg*t + fi)/(C*L*omg^2 - 1) + ... k1*sin(t/(sqrt(C)*sqrt(L)))

Obě analytická řešení si vzájemně odpovídají. Liší se pouze drobnými rozdíly v použitém značení. 47 Symbolem …“ je níže označeno pokračování kódu na dalším řádku. Do příkazové řádky se tečky nepíší a kód ” se nerozděluje.

86

7.2.3 Analytické řešení parciální diferenciální rovnice Pro úlohy definované na jednoduchých geometrických oblastech je možné s využitím integrálních transformací (Fourierova, Laplaceova) nalézt analytické řešení parciálních diferenciálních rovnic. Často uváděným příkladem postupu analytického řešení parciálních diferenciálních rovnic jsou rovnice vedení tepla v polonekonečné tyči nebo rovnice difúze. Tyto jevy jsou popsány parabolickou parciální diferenciální rovnicí, která je analogická například popisu skinefektu. Níže prezentovaný příklad popisuje vývoj výšky hladiny podzemní vody v blízkosti řeky, ve které vlivem přívalových dešťů krátkodobě stoupne hladina.

PŘÍKLAD 32 Hladina v řece Nise ovlivňuje výšku hladiny podzemní vody φ ve svém bezprostředním okolí. Tento vliv je možné matematicky popsat parabolickou parciální diferenciální rovnicí β

∂φ ∂2φ = ∂t ∂x2

s parametrem β, který reprezentuje difúzní materiálové vlastnosti horniny u břehů řeky. Ustálená, piezometrická výška hladiny podzemní vody u Nisy má hodnotu φ(x, 0) = 0 m. Uvedená hodnota představuje počáteční podmínku (chybí φ′ (x, 0)) příslušející výše zapsané rovnici. Vlivem přívalových dešťů se v čase t0 výška vodní hladiny v řece zvedne o 1 m. V čase t1 dojde k opětovnému poklesu hladiny v řece na původní úroveň. Výšky hladiny v řece definují jednu okrajovou podmínku φ(0, t) = H(t−t0 )−H(t−t1 ), kde H(. . .) je Heavisideova funkce a uvedený rozdíl představuje v čase obdélníkový puls. Druhou okrajovou podmínkou je hypotetická výška hladiny podzemní vody v nekonečnu φ(∞, t) = 0. Následuje souhrn informací známých ze zadání. 2

∂ φ β ∂φ ∂t = ∂x2 φ(x, 0) = 0 φ(0, t) = H(t − t0 ) − H(t − t1 ),

φ(∞, t) = 0

K řešení úlohy využijeme nejprve Laplaceovu transformaci. L [φ(x, t)] = ϕ(x, s) L [φ′ (x, t)] = s ϕ(x, s) − φ(x, [ ]0) s ϕ(x, s) − φ(x, 0) = β1 L

∂2φ ∂x2

=

1 ∂2ϕ β ∂x2

Pomocí transformace jsme získali obyčejnou diferenciální rovnici druhého řádu. Níže je z charakteristické rovnice odvozeno její obecné řešení. √ s = β1 λ2 , λ1,2 = ± βs √ √ ϕ(x, s) = c1 (s) e− βs x + c2 (s) e βs x Z okrajové podmínky φ(∞, t) = 0 plyne rovnost c2 (s) = 0. Zavedením substituce pro počáteční podmínku φ(x, 0) = φ0 , její následnou transformací ϕ(x, 0) = L [φ(x, 0)] =

φ0 = ϕ0 s

a zahrnutím do řešení uvažované diferenciální rovnice získáme předpis ve tvaru ϕ(x, s) = c1 (s) e−

87

√ βs x

+ ϕ0 .

Funkci c1 (s) vypočteme za pomoci okrajové podmínky φ(0, t). φ(0, t) = H(t − t0 ) − H(t − t1 ) −t s −t1 s ϕ(0, s) = e 0 −e s e−t0 s −e−t1 s = c1 (s) + ϕ0 s e−t0 s −e−t1 s − ϕ0 c1 (s) = s S vypočtenou funkcí c1 (s) už snadno získáme Laplaceův obraz ϕ(x, s) řešení zadané parabolické diferenciální rovnice. ( −t0 s ) √ e − e−t1 s ϕ(x, s) = − ϕ0 e− βs x + ϕ0 s Ke zpětné Laplaceově transformaci L−1 [ϕ(x, s)] = φ(x, t) jsou využity vzorce [ ] L−1 [ϕ0 ] = L−1 φs0 = φ0 , −cs L−1 [e [ −a√Fs ](s)] = H(t ( − c)f ) (t − c), ( ) a a e −1 √ √ L = erf c = 1 − erf , s 2 t 2 t kde erf (. . .) a erf c(. . .) značí po řadě chybovou funkci a doplňkovou chybovou funkci. Před využitím transformačních vzorců je vhodné upravit obraz řešení ϕ(x, s) do tvaru ϕ(x, s) = e

−t0 s e

√ − βs x

s

−e

−t1 s e

√ − βs x

s

− φ0

e−

√ βs x

s

+

φ0 s

Protože φ0 = 0, získáme zpětnou Laplaceovou transformací L−1 [ϕ(x, s)] = φ(x, t) funkci )) ( ( √ )) ( ( √ βx βx √ √ − H(t − t1 ) erf c , φ(x, t) = H(t − t0 ) erf c 2 |t − t0 | 2 |t − t1 | kde H(. . .) značí opět Heavisidovu funkci. Do výsledného vzorce bylo potřeba přidat absolutní hodnoty kvůli odmocninám ve jmenovatelích argumentů doplňkových chybových funkcí erf c(. . .). Výsledek odpovídá vzorcům uváděným v literatuře. Grafické znázornění výsledků ve vybraném čase bude kvůli srovnání provedeno v části 7.3.4 věnované numerickému řešení stejného příkladu.

7.3 Numerické řešení Analytické řešení matematických modelů je možné jen v jednoduchých případech a bývá uvedeno v rozsáhlejších monografiích. Často je nutno řešení získat numericky. Numerické řešení závisí na tom, jaký typ rovnic v modelu použijeme. Model lze popsat algebraickými, integrálními a diferenciálními rovnicemi nebo soustavami sestavenými z rovnic všech uvedených typů. Od druhu popisu se odvíjí způsob řešení. Následuje modelový příklad přechodového jevu v RL obvodu. Na tomto a na několika dalších příkladech bude dále v kapitole demonstrováno využití numerických metod řešení rovnic.

PŘÍKLAD 33 Zadání tohoto příkladu je převzato ze skript [16]. Cívka má indukčnost L a ohmický odpor RL . Cívka je zapojena do série s odporem R1 a paralelně k odporu R2 . Celá kombinace je připojena přes spínač na stejnosměrné napětí UV . Schéma obvodu ukazuje obrázek 38. Kotva relé přitáhne při proudu iL1 . Za jakou dobu dobu t1 relé po zapnutí kotvu přitáhne? 88

Parametry obvodu jsou: UV = 24 V, L = 6 H, RL = 200 Ω, R1 = 600 Ω, R2 = 300 Ω, iL1 = 18 mA.

Obrázek 38: Schéma RL-obvodu. Následuje přímé odvození funkčního předpisu časového vývoje proudu iL (t) procházejícího cívkou. Napětí zdroje se po sepnutí spínače rozloží na odporu R1 a na paralelním spojení cívky s odporem R2 . UV = UR1 + U2 = R1 i1 + i1 /G, Proud i1 tekoucí odporem R1 se ve větvích paralelního zapojení rozdělí na proudy i2 a iL v poměru dílčích vodivostí G2 , GRL . i1 = i2 + iL , G2 /GRL = i2 /iL , G = G2 + GRL i1 = iL (1 + G2 /GRL ) , UV = iL (1 + G2 /GRL ) (R1 + 1/(G2 + GRL )) Vodivosti G2 , GRL paralelně propojených částí obvodu je možné vyjádřit z následujících vzorců. U2 = UR2 = R2 i2 = i2 /G2 , U2 = UL = RL iL + L (diL /dt) = iL /GRL Vyjádřením vodivosti G = G2 + GRL a jejím dosazením do vztahu popisujícího součet napětí v obvodu získáme diferenciální rovnici s neznámou funkcí iL (t). UV

1 1 RL (R1 + R2 ) + R1 R2 R2 diL = iL + R1 + R2 L L R1 + R2 dt

Se zadanými parametry součástek v obvodu a s počátečními podmínkami lze nalézt řešení odvozené diferenciální rovnice ∂iL 4 200 iL (t) + = . (156) 3 ∂t 3 Řešení uvedené rovnice bylo provedeno v Sage. Kód následuje. sage: var('a b t') (a, b, t) sage: i_L = function('i_L',t) sage: deq = diff(i_L, t, 1) + a*i_L == b sage: proud = desolve(deq, i_L, [0,0], ivar=t) sage: var('U_V R_1 R_2 R_L L') (U_V, R_1, R_2, R_L, L) sage: prubeh = proud.substitute(a = R_L/L + R_1*R_2/(R_1 + R_2)/L, ...

89

b = U_V*R_2/(R_1 + R_2)/L); latex(prubeh) . . sage: U_V = 24 sage: L = 6 sage: R_L = 200 sage: R_1 = 600 sage: R_2 = 300 sage: prubeh = proud.substitute(a = R_L/L + R_1*R_2/(R_1 + R_2)/L, ... b = U_V*R_2/(R_1 + R_2)/L); latex(prubeh) . . sage: graf = plot(prubeh,0,0.1, color='red', title='Prechodovy jev', ... axes_labels=['$t[s]$','$i[A]$']) sage: graf.show() sage: graf.save('prechodovka.png')

Řešení sestavené diferenciální rovnice je v kódu uloženo do proměnné prubeh“. Výsledkem ” je funkce iL (t), která má tvar ) ) 200 200 1 ( ( 200 t) 1 ( iL (t) = e 3 − 1 e(− 3 t) = 1 − e(− 3 t) . 50 50 Odtud už se dá vypočítat čas kdy dosáhne proud tekoucí cívkou hodnoty 18 mA nutné pro přitažení kotvy. t = 0,0345387763949107 s ∼ = 35 ms Výsledek lze orientačně ověřit z grafu průběhu funkce iL (t) v čase na obrázku 39.

Obrázek 39: Časový vývoj proudu tekoucího cívkou. Vypočítaná hodnota času přitažení kotvy je řešením nelineární rovnice ) 200 1 ( 1 − e( − 3 t ) 0,018 = 50 90

a poslouží jako referenční řešení pro srovnání výsledků níže představených numerických metod. Diferenciální rovnice (156) s hledanou, neznámou funkcí iL (t) bude dále použita také k předvedení způsobu použití metod numerické integrace. 7.3.1 Metody numerického řešení algebraických a transcendentních rovnic Pokud nemá matematický popis vybraného jevu speciální, jednoduchý tvar, potom není snadné a často ani možné najít analytické řešení použitých rovnic. To platí například už pro formálně jednoduché algebraické rovnice s polynomy vyšších řádů. V této části textu budou představeny dva jednoduché postupy řešení algebraických a transcendentních rovnic. Budou to bisekce (metoda půlení intervalu) a Newtonova metoda. Oba postupy budou aplikovány na řešení rovnice 0=

) 200 1 ( 1 − e(− 3 t) − 0,018 50

popisující okamžik přitažení kotvy během přechodového jevu v RL obvodu (viz výše uvedený příklad). Metoda půlení intervalu - bisekce

PŘÍKLAD 34 K vysvětlení principu metody půlení intervalu poslouží grafická interpretace významu rovnice. Mějme dánu dříve připravenou rovnici 0=

) 200 1 ( 1 − e(− 3 t) − 0,018. 50

Na její pravé straně stojí funkční předpis f (t) a na straně levé potom jedna speciální hodnota funkce, hodnota 0. Kořen rovnice je místem, kde funkce f (t) protíná osu s nezávisle proměnnou, v tomto případě t (viz obrázek 40). Z grafu je zřejmé, že hledaný průsečík leží v intervalu t ∈ ⟨0,02; 0,04⟩. Meze intervalu byly v tomto případě zvoleny tak, aby funkční hodnoty v krajních bodech měly různá znaménka. To je u spojité funkce f (t) dostačující podmínka k tomu, aby vybraný interval obsahoval její kořen. Pokud jsme tázáni na řešení rovnice, pak zpravidla nepostačí odpověď ve formě určení intervalu, jehož meze se neshodují ani v první platné číslici. Když chceme odhad řešení zlepšit, musíme interval zmenšovat do okamžiku než dosáhneme požadované přesnosti řešení.

91

Obrázek 40: Hledaný průsečík funkce f (t) s osou t. Metoda půlení intervalu je postupem systematického přibližování mezí intervalu obsahujícího kořen řešené rovnice. Její postup je následující: 1. Určíme počáteční odhad mezí intervalu řešení ⟨a; b⟩. V našem případě byl zvolen interval t ∈ ⟨0,02; 0,04⟩. V krajních bodech intevalu musí mít funkce f (t) různá znaménka. 2. Vyhodnotíme znaménko funkce v polovině zvoleného intervalu, tj. f ((b + a)/2). V příkladu nás zajímá znaménko hodnoty funkce f (t) v bodě t = 0,03. f (0,03) < 0 3. Vybereme a nahradíme tu z mezí intervalu ve které má funkce f (t) stejné znaménko jako v bodě f ((b + a)/2). Interval se tím zmenší na polovinu. Hodnota a nebo hodnota b se změní. Pro příklad s kotvou to znamená t ∈ ⟨0,03; 0,04⟩. 4. Dále určíme znaménko funkce f (t) v polovině nového intervalu ⟨a; b⟩. Pro příklad RL obvodu vyjde f (0,35) > 0, takže nahradíme horní mez intervalu a dostaneme t ∈ ⟨0,3; 0,35⟩. 5. Zjistíme znaménko funkce f (t) v polovině nového intervalu ⟨a; b⟩. Protože f (0,325) < 0, tak nahradíme dolní mez a dále budeme hledat řešení uvnitř intervalu ⟨0,325; 0,35⟩. 6. . . . 7. S půlením pokračujeme do chvíle, kdy vzdálenost mezí intervalu klesne pod dvojnásobek předem zadané přesnosti řešení |b − a| <= 2 · ε.

92

a 0,020000 0,030000 0,030000 0,032500 0,033750 0,034375 0,034375 0,03453125 0,03453125 0,03453125 0,03453125

b 0,040000 0,040000 0,035000 0,035000 0,035000 0,035000 0,03468750000 0,03468750000 0,03460937500 0,03457031250 0,03455078125

(b + a)/2 0,030000 0,035000 0,032500 0,033750 0,034375 0,03468750000 0,03453125000 0,03460937500 0,03457031250 0,03455078125 0,034541015625

znaménko f ((b + a)/2) − + − − − + − + + +

Numericky získaný výsledek má hodnotu t = 3,4541015625 · 10−2 ± 9,765624 · 10−6 s. Pro srovnání připomeňme výsledek analytického výpočtu stejné úlohy, který byl t = 0, 0345387763949107 s. S ohledem na zadanou přesnost výpočtu jsme dosáhli dobré shody výsledků. Doposud bylo při ukázkách řešení modelových příkladů předpokládáno zadávání příkazů do konzole (příkazové řádky) programu Sage. U delšího, strukturovaného kódu je vhodnějším řešením příprava skriptu s příponou .sage“. Níže uvedené příkazy byly zapsány do skriptu ” bisekce.sage“ a volány příkazem load(’./bisekce.sage’)“. ” ” var('t') prechod = 1/50 * (1 - e^(-200*t/3)) - 0.018 hodnota = fast_float(prechod, 't') def bisekce(f, a, b, chyba): while (b-a).abs()>=(2*chyba): c = (a + b)/2 if f(c) < 0.0: a = c if f(c) > 0.0: b = c if f(c) == 0.0: return c print a, b, (b+a)/2, (b-a)/2 return (b+a)/2 bisekce(hodnota, 0.02, 0.04, 0.00001)

Newtonova metoda (metoda tečen) Jiným postupem pro řešení algebraických a transcendentních rovnic je iterační Newtonova metoda (metoda tečen).

PŘÍKLAD 35 Uvažujme stejný scénář jako v příkladu řešeném metodou půlení intervalu. ( ) 200 1 iL (t) = 50 1 − e( − 3 t ) ( ) 200 f (t) = 1 1 − e(− 3 t) − 0,018 = 0,0 50

RL obvodem protéká proud. Zajímá nás čas od sepnutí spínače kdy proud tekoucí cívkou dosáhne velikosti 18 mA a vyvolá přitažení kotvy relé. 93

Pro řešení rovnice f (t) = 0 Newtonovou metodou potřebujeme znát (určit) derivaci f ′ (t) a navíc musíme mít (udělat) počáteční odhad řešení t0 . 4 (− 200 t) e 3 3 Znalost derivace f ′ (t) a počátečního odhadu řešení t0 využijeme k postupnému zpřesňování hodnoty kořene řešené rovnice. Zpřesňování probíhá podle schématu daného vzorcem f ′ (t) =

tk+1 = tk −

f (tk ) f ′ (tk )

Uvedený vzorec říká: Vydej se z bodu [tk , f (tk )] ve směru tečny k ose t. Místo, kde se zkříží ” tvoje cesta s osou t je hledané tk+1 .“ :) Zpřesňování odhadu řešení se provádí do chvíle než velikost funkční hodnoty f (t) v místě aktuálního odhadu řešení tk klesne pod stanovenou mez přesnosti ε, tj. |f (tk )| < ε. Za pomoci představeného postupu byla při zadaném počátečním odhadu řešení t0 = 0,025 s daná rovnice numericky vyřešena ve třech krocích. cas tk 0,0320582649242949 0,0343445316457873 0,0345375241055565

funkcni hodnota f (tk ) −0, 000359653088736 −2,60677195233 · 10−5 −1,66978883988 · 10−7

derivace f ′ (tk ) 0,157310205916 0,135071181302 0,133344465259

Nalezený numerický odhad řešení t = 0,0345375241055565 s je blízký analytickému řešení t = 0,0345387763949107 s. Rychlost (resp. počet kroků zpřesnění) odhadu řešení rovnice Newtonovou metodou závisí na stanovené přesnosti řešení a na vzdálenosti počátečního odhadu řešení od jeho skutečné hodnoty. Závislost délky (doby trvání) výpočtu na volbě počátečního odhadu řešení t0 si mohou zájemci ověřit spuštěním níže uvedeného skriptu. var('t') prechod = 1/50 * (1 - e^(-200*t/3)) - 0.018 def newtonova_metoda(f, vychozi_bod, tolerance): hodnota_funkce = fast_float(prechod, 't') derivace = diff(f, t, 1) hodnota_derivace = fast_float(derivace, 't') while f(vychozi_bod).abs() > tolerance: vychozi_bod -= hodnota_funkce(vychozi_bod)/hodnota_derivace(vychozi_bod) print vychozi_bod, hodnota_funkce(vychozi_bod), hodnota_derivace(vychozi_bod) return vychozi_bod print "Vychozi bod t = 0.025" newtonova_metoda(prechod, 0.025, 0.00001) print "Vychozi bod t = 0.1" newtonova_metoda(prechod, 0.1, 0.00001)

Při použití Newtonovy metody se mohou objevit komplikace v případech, kdy f ′ (xk ) = 0 nebo když má funkce lokální extrém a odhad řešení kolem extrému osciluje mezi dvěma body, tj. pokud f (tk+1 ) f (tk ) tk+1 − ′ = tk − ′ . f (tk+1 ) f (tk ) V těchto situacích nevede aplikace Newtonovy metody k nalezení řešení rovnice. 94

7.3.2 Základní iterační metody řešení soustav lineárních algebraických rovnic Principy numerického řešení soustav lineárních algebraických rovnic (SLAR) budou níže ukázány na řešení zvolené modelové úlohy dvěma základními metodami. Příklad ilustruje skutečnost, že využití numerických metod nemusí vždy vést k získání hledaného řešení. Aby se numerická aproximace blížila řešení soustavy, musí soustava zpravidla splňovat určitá kritéria. Známými metodami řešení soustav lineárních algebraických rovnic jsou Jacobiho a GaussSeidlova metoda. Obě metody vycházejí ze speciálního rozdělení matice A na diagonální (D), horní trojúhelníkovou (U) a dolní trojúhelníkovou (L) část. A = D + L + U,

→ − → − − → A− x = D→ x + L→ x + U− x = b

Vzorce pro řešení soustav s využitím vybraných metod jsou uvedeny v následujícím výčtu. • Jacobiho metoda (Metoda prosté iterace) (− )) ) →(k) → ( − → ( − → − → → → x (k+1) = D−1 b − (L + U) − x (k) x (k) = D−1 b − D−1 (L + U) − x = D−1 b − T −

− − má vyšší paměťové nároky. Je třeba ukládat samostatně → x (k) i → x (k+1) . Nedochází k prů→ běžné aktualizaci vektoru − x (k) . • Gaus-Seidlova metoda

(− → ( − (k) )) − → − → → x (k+1) = (D + L)−1 b − U → x = (D + L)−1 b − T − x (k) ,

aktualizuje složky stavového vektoru průběžně. Postačující podmínkou konvergence výše uvedených metod je spektrální poloměr (v absolutní hodnotě největší vlastní číslo ρ(T) = max{|λ1 |, . . . , |λn |}) iterační matice T menší než 1, zapsáno algebraicky ρ(T) < 1.

PŘÍKLAD 36 Připomeňme si příklad výpočtu proudů v obvodu tvořeném dvěma zdroji stejnosměrného napětí a třemi rezistory (viz obrázek 41). V jedné z předchozích kapitol byl na tomto příkladu předveden způsob symbolického řešení soustavy rovnic v programu Sage.

Obrázek 41: Modelový obvod. Soustava pro výpočet proudů (I1 , I2 , I3 ) tekoucích obvodem má tvar −I1 + I2 − I3 = 0, R1 I1 + R2 I2 = U1 , R2 I2 + R3 I3 = U2 .

95

Rozšířený maticový zápis soustavy rovnic s dosazenými parametry R1 = 50 Ω, R2 = 75 Ω, R3 = 100 Ω, U1 = 12 V, U2 = 24 V, má tvar



 −1 1 −1 0  50 75 0 12  . 0 75 100 24

Pro iterační matici Jacobiho metody je podmínka konvergence splněna (viz výstupy níže uvedeného kódu). Iterační matice Gauss-Seidlovy metody podmínku konvergence nesplňuje. Následuje kód skriptu pro program Sage. reset() print 'Kirchhoffovy zakony - Jacobiho metoda' R_1 = 50.0 R_2 = 75.0 R_3 = 100.0 U_1 = 12.0 U_2 = 24.0 A = Matrix([[-1, 1, -1],[R_1, R_2, 0],[0, R_2, R_3]]); A b = vector([0, U_1, U_2]); b D = diagonal_matrix(A.diagonal()); D LU = A-D; LU D_inverze = D.inverse(); D_inverze n_iteraci = 500; odhad = vector([0.2, 0.2, 0.2]) odhad2 = odhad while n_iteraci > 1: odhad2 = D_inverze*(b - LU*odhad) n_iteraci -= 1 odhad = odhad2 T1 = matrix(QQ, D.inverse()*(LU)) print 'Vlastni cisla' print T1.eigenvalues() print 'Reseni' print A.inverse()*b print odhad reset() print 'Kirchhoffovy zakony - Gauss-Seidlova metoda' import numpy R_1 = 50.0 R_2 = 75.0 R_3 = 100.0 U_1 = 12.0 U_2 = 24.0 A = Matrix([[-1, 1, -1],[R_1, R_2, 0],[0, R_2, R_3]]) b = vector([0, U_1, U_2]) U = Matrix(numpy.triu(A)) - diagonal_matrix(A.diagonal()) DL = Matrix(A-U)

96

DL_inverze = DL.inverse() n_iteraci = 500; odhad = vector([0.2, 0.2, 0.2]) while n_iteraci > 1: odhad = DL_inverze*(b - U*odhad) n_iteraci -= 1 T2 = matrix(QQ, DL.inverse()*U) print 'Vlastni cisla' print T2.eigenvalues() print 'Reseni' print A.inverse()*b print odhad

Díky splnění postačující podmínky konvergoval odhad řešení Jacobiho metodou [I1 ; I2 ; I3 ] = [0,0184616363349266; 0,147692326275641; 0,129230709714611] A k řešení analytickému. [I1 ; I2 ; I3 ] = [0,0184615384615385; 0,147692307692308; 0,129230769230769] A U Gauss-Seidlovy metody ke konvergenci numerického odhadu k analytickému řešení nedošlo. Výsledkem numerického odhadu byly proudy [I1 ; I2 ; I3 ] = [−4,03430868754227 · 1031 ; 2,68953912502818 · 1031 ; −2,01715434377114 · 1031 ] A

7.3.3 Numerické řešení obyčejných diferenciálních rovnic V praxi mohou nastat případy, kdy jsme postaveni před úkol najít primitivní funkci (neurčitý integrál), která není známa. Znamená to, že nedokážeme získat funkční předpis integrálu, ale nic není ztraceno. V takovém případě zbývá možnost odhadnout určitý integrál funkce na zadaném intervalu, tj. přibližně určit plochu, která je v určeném intervalu pod křivkou funkce (pro jednoduchost uvažujeme případ, kdy funkce nenabývá záporných hodnot, tj. f (. . .) ≥ 0).

PŘÍKLAD 37 Výpočet obsahu plochy pod křivkou funkce je níže proveden několika metodami na už opakovaně použitém příkladu RL obvodu (viz obr. 38). Obrázek 42 ukazuje graf funkce, která je vyjadřením vývoje časové derivace proudu tekoucího cívkou v uvažovaném obvodu.

97

Obrázek 42: Graf časového vývoje derivace proudu tekoucího cívkou zařazenou v RL obvodu. Diferenciální rovnice

∂iL (t) = 1,33 ∂t popisující proud tekoucí cívkou je níže řešena Eulerovou metodou a metodou Runge-Kutta. Kromě toho je proveden dolní, horní a střední odhad plochy pod křivkou pomocí Obdélníkového ” pravidla“. Všechna řešení jsou realizovatelná spuštěním níže vloženého kódu pro prostředí Sage. Tvar zadání diferenciální rovnice se v Sage mírně liší podle zvolené metody řešení. 66,67 iL (t) +

reset() var('a b t') i_L = function('i_L',t) deq = diff(i_L, t, 1) + a*i_L == b U_V = 24.0 L = 6.0 R_L = 200.0 R_1 = 600.0 R_2 = 300.0 deq_s_parametry = deq.substitute(a = R_L/L + R_1*R_2/(R_1 + R_2)/L, ... b = U_V*R_2/(R_1 + R_2)/L) print deq_s_parametry # print 'Analyticke reseni' # print '-----------------' reseni_ana = desolve(deq_s_parametry, i_L, [0.0,0.0], ivar=t) # print 'Eulerova metoda' # print '---------------' from sage.calculus.desolvers import eulers_method t_eu,i_Leu = PolynomialRing(RR,2,"i_L").gens() deq_eu = U_V*R_2/(R_1 + R_2)/L - (R_L/L + R_1*R_2/(R_1 + R_2)/L)*i_Leu

98

reseni_eu = eulers_method(deq_eu, 0, 0, 0.02, 0.1, algorithm='None') # print 'Metoda Runge-Kutta 4. radu' # print '--------------------------' from sage.calculus.desolvers import desolve_rk4 reseni_RK4 = desolve_rk4(deq_s_parametry, i_L, ics=[0,0], end_points=0.1, step=0.02) print 'Plocha pod krivkou zadane funkce' print '---------------------------------------' der_i_L(t) = b*exp((-1)*a*t) reseni_ui = (integral(der_i_L, t, 0, 0.1)).substitute(a = R_L/L + ... R_1*R_2/(R_1 + R_2)/L, b = U_V*R_2/(R_1 + R_2)/L) hodnota_der = fast_float((der_i_L.substitute(a = R_L/L + ... R_1*R_2/(R_1 + R_2)/L, b = U_V*R_2/(R_1 + R_2)/L)), 't') def urcite_integraly(funkce, zacatek, konec, pocet_podintervalu): dolni_odhad = 0.0 prostredni_odhad = 0.0 horni_odhad = 0.0 pozice = zacatek krok = (konec-zacatek)/pocet_podintervalu while pozice < (konec-krok): dolni_odhad += krok*(funkce(pozice)) prostredni_odhad += krok*(funkce(pozice+krok/2)) horni_odhad += krok*(funkce(pozice+krok)) pozice += krok; return [dolni_odhad, prostredni_odhad, horni_odhad] [dolni, stredni, horni] = urcite_integraly(hodnota_der, 0.0, 0.1, 50) print print print print print print

'dolni odhad urciteho integralu' dolni 'prostredni odhad urciteho integralu' stredni 'horni odhad urciteho integralu' horni

# print 'Srovnani vysledku' # print '-----------------' graf_srovnani = plot(reseni_ana, 0, 0.1, color='blue', ... title='Srovnani numerickych metod', axes_labels=['$t[s]$','$i[A]$'], ... legend_label='analyticky') + list_plot(reseni_eu, color='red', ... plotjoined=True, legend_label='eulerova_metoda') + ... list_plot(reseni_RK4, color='green', plotjoined=True, legend_label='metoda RK4') graf_srovnani.show() graf_srovnani.save('./srovnani_metod.png')

Výsledkem běhu uvedeného skriptu je graf na obrázku 43. Porovnání výsledků ukazuje, že metoda vyššího řádu (Runge-Kutta 4. řádu) dosahuje při délce časového kroku 0,02 s lepších výsledků než Eulerova metoda. Obrázek 44 ukazuje zlepšení výsledků obou metod dosažené zmenšením délky časového kroku na polovinu původního, tj. na 0,01 s.

99

Obrázek 43: Výsledky numerické integrace s časovým krokem 0,02 s.

Obrázek 44: Výsledky numerické integrace s časovým krokem 0,01 s.

100

7.3.4 Numerické řešení parciálních diferenciálních rovnic Kvůli obtížnosti analytického řešení parciálních diferenciálních rovnic popisujících reálné jevy jsou k odhadu jejich řešení často využívány numerické metody. V technické praxi nalézá dnes široké uplatnění metoda konečných prvků, ale vzhledem k větší názornosti bude níže k řešení modelového příkladu použita metoda konečných diferencí (metoda sítí). Metoda konečných diferencí nahrazuje derivace rozdíly hodnot hledané veličiny v uzlových bodech diskretizační sítě. Tyto rozdíly jsou váženy délkou prostorového nebo časového kroku, podle druhu nezávisle proměnné. Diskretizace úlohy v prostoru vede na soustavu obyčejných diferenciálních rovnic, které nahradí původní parciální diferenciální rovnici.

PŘÍKLAD 38 V kapitole věnované analytickému řešení parciálních diferenciálních rovnic byl uveden modelový příklad řešení parabolické, parciální diferenciální rovnice popisující vývoj výšky hladiny podzemní vody (φ(x, t)) v okolí řeky postižené nárazovou povodní. Pro připomenutí následuje matematický zápis zadání úlohy. 2

∂ φ β ∂φ ∂t = ∂x2 φ(x, 0) = 0 φ(0, t) = H(t − t0 ) − H(t − t1 ),

φ(∞, t) = 0

K numerickému řešení je třeba diskretizovat rovnici v prostoru a v čase. Simulovaný časový interval 2,0 hodiny byl v řešené úloze rozdělen na nout = 2000 časových kroků. Délka trvání povodně byla doba_povodne = 0,3 hodiny. Pro porovnání výsledků analytického a numerického výpočtu byl zvolen čas vyber_cas = 1,75 hodiny. Sledovaná jednorozměrná oblast dlouhá xmax = 5 m byla rozdělena na 30 prostorových kroků (29 vnitřních uzlů). K prostorové diskretizaci úlohy posloužilo centrální schéma (viz níže vložený kód skriptu pro Sage). K vyřešení v prostoru diskretizované úlohy posloužila implementace metody Runge-Kutta 4. řádu obsažená v Sage. Kód příslušného skriptu následuje. reset() from sage.calculus.desolvers import desolve_system_rk4 x,t = var('x t') xmax = 5 # metru mx = 29 # 29 vnitrnich uzlu, tj. 30 dilu krok = xmax/(mx + 1) t0 = 0.0 # zacatek simulovaneho intervalu konec = 2.0 # 2.0 hodiny - konec simulovaneho intervalu nout = 2000 # pocet casovych kroku simulace, vcetne pocatecni podminky dt = konec/nout # delka casoveho kroku doba_povodne = 0.3 # hodiny, tj. 18 minut, kratka povoden vyber_cas = 1.75 # 1hod 45 minut cas vyhodnocovani rozlozeni vysek hladiny vc_idx = round(vyber_cas/dt) if vc_idx > (nout - 1): print 'vybrany cas nelezi v simulovanem intervalu' print vyber_cas, vc_idx, konec betka2 = 1e+2 # parametr difuze, (storativita/transmisivita)^2, (1e-4/1e-5)^2 ict = 0 # piezometricka vyska v metrech bcd = 1*(sign(t-t0)/2 - sign(t-t0-doba_povodne)/2) # pohyb hladiny (obdelnikovy puls) # definice matice promennych hodnot hladin v uzlech

101

hladiny = [var('h_%d' % i) for i in range(mx)] pocatecni_podminky = [ict for i in range(mx+1)] pocatecni_podminky[0] = t0; # prvni hodnota je cas pocatku simulovaneho intervalu, t0 # definovani soustavy ODE vznikle diskretizaci PDE v prostoru de = [t for i in range(mx)] # jde o to aby byly definovany funkce casu for idx in range(mx): if idx == 0: de[idx] = (1/betka2)*(hladiny[idx+1] - 2*hladiny[idx] + bcd)/krok^2 elif idx == (mx - 1): de[idx] = (1/betka2)*(hladiny[idx] - 2*hladiny[idx] + hladiny[idx-1])/krok^2 else: de[idx] = (1/betka2)*(hladiny[idx+1] - 2*hladiny[idx] + hladiny[idx-1])/krok^2 # reseni, alternativne by slo pouzit desolve_system, nebo desolve_odeint para_reseni = desolve_system_rk4(de, hladiny, ics=pocatecni_podminky, ivar=t, ... end_points=konec, step=dt) # analyticke reseni pro porovnani analyt_res = heaviside(t-t0)*(1-erf(x/2/sqrt((1/betka2)*abs(t-t0)))) - ... heaviside(t-t0-doba_povodne)*(1-erf(x/2/sqrt((1/betka2)*abs(t-t0-doba_povodne)))) + ... ict #ROZLOZENI HLADIN V PROSTORU VE ZVOLENEM CASE analyt_rozlozeni = analyt_res.substitute(t=vyber_cas) # zobrazit prostorove rozlozeni vysky hladin ve zvolenem case vzdalenost = [i*krok for i in range(mx+1)] # okrajova podminka ve zvolenem case okrajova = bcd.substitute(t = vyber_cas) # prikaz zip mi umozni kombinovat pri tisku dva seznamy graf = list_plot(zip(vzdalenost, [okrajova] + para_reseni[vc_idx][1:]), ... plotjoined=True, color='red', title='Vysky hladiny', ... axes_labels=['$x[km]$','$h[m]$'], legend_label="numericky") graf2 = plot(analyt_rozlozeni, (x,0,xmax), color='blue', legend_label='analyticky') (graf + graf2).show() (graf + graf2).save('vyvoj_v_bode_srovnani.png')

Oba způsoby řešení, analytický i numerický přístup, vedly ke srovnatelným výsledkům. Graf na obrázku 45 ukazuje prostorové rozložení výšky hladiny podzemní vody v čase 1,75 hodiny od začátku simulovaného časového intervalu.

102

Obrázek 45: Prostorové rozložení výšky hladiny podzemní vody. V případě zjemňování prostorové diskretizace úlohy se výsledky analytického a numerického řešení modelu vzájemně přibližují. Cenou za zvýšení přesnosti výsledků numerického výpočtu je prodloužení výpočetního času modelu.

Klasifikace a význam okrajových podmínek Je známo, že úlohy popsané diferenciálními rovnicemi (obyčejnými i parciálními) se skládají jednak z příslušné rovnice, jednak z dodatečných podmínek, podle kontextu nazvaných okrajové nebo počáteční. Z přirozenosti věci plyne, že při zkoumání jevů v systému se neobejdeme bez specifikace interakce systému s okolím. V našem případě diferenciální rovnice popisuje fyzikální podstatu a okrajová podmínka interakci s okolím. Pro eliptické rovnice 2. řádu rozlišujeme tyto tři typy okrajových podmínek: • 1. druhu (Dirichletova) – předepsán je průběh potenciálu u = uD .

(157)

Pod symbolem uD může být např. teplota, elektrický potenciál, elastické posunutí apod. na zvolené ploše. Obecně to je funkce. Nejjednodušším a v praxi často používaným případem je konstantní hodnota, např. potenciál na elektrodě. Pokud je uzemněná, je nulový. • 2. druhu (Neumannova) – předepsán je průběh derivace potenciálu. − → → → (K ∇u) · − n = Kgrad(u) · − n = qN ,

(158)

→ kde − n je normála k okrajové ploše a K je koeficient. Fyzikální význam závisí na typu řešeného problému. U transportních jevů jde většinou o tok, např. tok tepelné energie z plošného zdroje tepla. V problémech elasticity to může být předepsaný průběh deformace či elastického napětí. 103

• 3. druhu (Cauchyova, Newtonova) – kombinace potenciálu a jeho derivace − → → → (K ∇u) · − n + λu = K · grad(u) · − n + λu = q3 .

(159)

Podmínku 3. druhu je nejpřirozenější chápat jako závislost toku přes hranici na rozdílu hodnoty potenciálu uvnitř (neznámá) a vně (referenční zadaná hodnota). Vztah lze ekvivalentně zapsat → − → → (K ∇u) · − n = K · grad(u) · − n = λ(u − u3 ), q3 = λu3 , (160) což je obvyklá forma pro zadání parametrů ve výpočetních softwarech. Pro hodnoty λ → 0 se podmínka transformuje na 2. druh, pro hodnoty λ → ∞ se podmínka transformuje na 1. druh, je tedy zobecněním výše uvedených. Podmínku lze s výhodou užít tam, kde bychom předepsali podmínku prvního typu, ale spíše ve smyslu doporučení“ než nařízení (přesného zadání)“ ” ” a hodnota λ vyjadřuje míru (váhu) určitosti zadání. V úlohách vedení tepla popisuje podmínka 3. druhu přestup tepla, poněvadž tepelný tok je úměrný rozdílu teploty povrchu tělesa a okolí. Tím se zahrnují jevy v tenké přípovrchové vrstvě. V úlohách mechaniky popisuje pružné uložení. Podmínky existence a jednoznačnosti řešení okrajových úloh Při konstrukci nějakého modelu v podobě okrajové úlohy pro diferenciální rovnici má smysl se ptát, zda takto zadaná úloha vůbec má nějaké řešení a zda toto řešení je jednoznačné, tj. zda nemůže existovat více různých řešení. Přirozená intuice nám sice říká, že když modelujeme nějaký reálný systém, tak ten přece existuje a má určité vlastnosti a ne jiné a to je tedy ono jednoznačné řešení naší úlohy. Problém je v tom, že na papíře nebo v počítači řešíme úlohu, která je jen modelem reality a takto zjednodušená úloha již nemusí mít všechny potřebné vlastnosti původního reálného systému (možná jsme zanedbali nějakou vlastnost, která je ve skutečnosti pro probíhající děje, tedy pro existenci správného řešení, klíčová). Z pohledu přesné matematické formulace je problematika existence a jednoznačnosti řešení docela složitá. Zde zmíníme jen hlavní aspekty, důležité při zadávání úloh do výpočetního software. Například přesné splnění diferenciálních rovnic a okrajových podmínek vyžaduje hodně přísné podmínky na spojitost neznámých funkcí (pole počítané veličiny), parametrů (rozložení materiálových vlastností) a tvaru oblasti, které reálně ani nemohou být splněny (těleso má rohy a hrany, na sebe navazují dva různé materiály apod.), a nelze je tedy požadovat při zadávání úloh do výpočetních programů. Některé matematické podmínky tedy nejsou tak podstatné pro získání rozumného řešení pomocí modelovacího software, některé však je nutno dodržet, nebo alespoň vědět, jaké případné závady v řešení při jejich nedodržení hrozí. Základní pravidlo pro jednoznačné řešení: u eliptické rovnice 2. řádu (tj. modely stacionárního vedení tepla, elektrostatiky, elasticity, atd.) musí být alespoň na jedné části hranice předepsána okrajová podmínka prvního nebo třetího druhu (tj. nelze formulovat úlohu jen s podmínkami druhého druhu). U úlohy s předepsanou podmínkou 2. druhu v ustáleném stavu nastávají dva případy: • Zadané toky dovnitř a ven z oblasti jsou v rovnováze“ se zdroji, tj.48 ”∫ ∫ qN dS − f dV = 0, Γ

Ω

potom má úloha nekonečně mnoho řešení, které se navzájem liší o aditivní konstantu. V úloze jsme všude zadali jen derivace, takže po přičtení konstanty zůstávají hodnoty všech derivací stejné. Hodnoty potenciálu tedy nejsou jednoznačně určeny49 , ale hodnoty toku (derivace) již jednoznačně určeny jsou. 48 49

Symbolem Ω je označen vyšetřovaný objem, Γ je jeho povrch (hranice). To je ovšem jeho obecná vlastnost, jak jsme již několikrát zdůraznili.

104

• Zadané toky dovnitř a ven z oblasti nejsou v rovnováze“ se zdroji, tj. ” ∫ ∫ qN dS − f dV ̸= 0, Ω

Γ

potom úloha nemá žádné řešení a nemá ani rozumný fyzikální smysl. V ustáleném stavu nesplňuje zákon zachování energie. Takové zadání nemohlo vzniknout přílišným zjednodušením reality, ale spíše zcela chybnou úvahou. Úlohy s předepsanou podmínkou 2. druhu splňující podmínku rovnováhy toků (podmínka kompatibility) tedy svůj význam mít mohou: např. známe celkový tok, zajímá nás jeho rozložení v prostoru, a přitom nás nezajímají hodnoty potenciálu. To, jak na takové zadání reaguje použitý software je různé: • Buď ohlásí, že je zadáno málo okrajových podmínek (nedourčená úloha) a odmítne výpočet provést, • nebo v tichosti spočítá výsledek s tím, že si automaticky použil nějakou předdefinovanou hodnotu potenciálu. Každopádně je však lepší mít jistotu co počítáme a tuto hodnotu si raději zvolit sami, tj. předepsat okrajovou podmínku 1. druhu v jednom zvoleném referenčním bodě. V úlohách mechaniky (pružnosti) odpovídá splnění uvedených podmínek tomu, že těleso je ve statické rovnováze (nemůže se vlivem zadaných sil začít pohybovat). Doplnění okrajové podmínky 1.druhu odpovídá přidání podpěry, v níž (jako výsledek řešení úlohy) působí nulová síla. V úlohách elektromagnetizmu vzniká další závažný problém s okrajovými podmínkami v nekonečnu. Zde musí být vypočtené hodnoty nulové. Nejsnáze se obchází tím, že počítaná oblast se zvětší tak, aby na její hranici byla tato podmínka přibližně splněna. Pak se na její hranici dosadí nulová hodnota. Prakticky to znamená, že hranice je ve vzdálenosti asi desetinásobné vzhledem k nejdelšímu geometrickému rozměru modelu. Rozměry studované oblasti se tedy výrazně (a zbytečně) zvětší.

7.4 Metody zpracování naměřených dat Někdy je hledaným modelem funkční předpis mezi jednotlivými sledovanými veličinami. Takový → − funkční předpis lze v případě existující závislosti získat z naměřených dat {− x,→ y }. − → x − → y

x0 y0

x1 y1

... ...

xk−1 yk−1

Mohou nám k tomu pomoci následující postupy: • Interpolace – Lagrangeova interpolace. – Splinovací funkce. • Regrese - metoda nejmenších čtverců. 7.4.1 Interpolace V případě interpolace naměřených dat neuvažujeme existenci chyb měření. Výsledná funkce prochází měřením definovanými body.

105

Lagrangeova interpolace Provádíme-li Lagrangeovu interpolaci, pak dochází k prokládání všech bodů jedním polynomem. Lagrangeův interpolační polynom vycházející ze znalosti n měření má stupeň n − 1. Pomocí sumace je soustava rovnic pro konstrukci Lagrangeova interpolačního polynomu N (x) vyjádřena vzorci j−1 n ∑ ∏ N (x) = aj lj (x) − yj , lj (x) = (x − xi ), j=0

i=0

kde neznámé jsou koeficienty aj . Proměnné xj , yj reprezentují naměřené hodnoty. V maticové formě jde soustavu pro výpočet koeficintů aj formujících Lagrangeův interpolační polynom v Newtonově tvaru zapsat výrazně srozumitelněji v následující podobě       1 0 ... ... ... y0 a0  1 (x − x0 )  0 . . . . . .    a 1   y1    1 (x − x0 ) (x − x0 )(x − x1 ) 0     . . .    ..   ..        .. .. ..  1 · . = .  . . . 0    ..   ..  k−1    .   .  ∏   1 (x − x0 ) (x − x0 )(x − x1 ) . . . (x − xi ) yk−1 ak−1 i=0

Výhodou uvedeného zápisu je možnost přidávat body bez nutnosti kompletního přepočítávání koeficientů aj . Dále je na příkladu demonstrována zavedená metoda výpočtu koeficientů. Postup: 1. Naměřená data uspořádáme do tabulky a řádkům tabulky přiřadíme indexy i ∈ {0, . . . , n − 1}. 2. Vypočteme postupně diference f (x0 ), . . . , f (x0 , x1 , . . . , xn−1 ) nultého až (n − 1)-ního řádu. Diference nultého řádu je samotná hodnota funkce v bodě f (xi ), i ∈ {0, . . . , n − 1}. Diference řádu k se počítá z diferencí nižšího řádu podle vzorce f(xi , . . . , xi+k+1 ) =

f(xi+1 , . . . , xi+k+1 ) − f(xi , . . . , xi+k ) xi+k+1 − xi

3. Počítané diference je zvykem průběžně psát do nových sloupců dříve vytvořené tabulky vždy od k-tého řádku dolů. 4. Sestavíme polynom N (x) = f0 (x0 ) + f1 (x0 , x1 )(x − x0 ) + f2 (x0 , x1 , x2 )(x − x0 )(x − x1 ) + . . . . . . + fn−1 (x0 , . . . , xn−1 )(x − x0 ) · . . . · (x − xn−2 ), kde fj (. . .), j ∈ {0, . . . , n − 1} stojící v tabulce na diagonále (tj. na vrcholu sloupců) je hledané aj .

PŘÍKLAD 39 Zadání a část řešení byly převzato z [14]. Jsou dány body [−2, −9][0, 1][1, 3][2, 11]. Určerte Lagrangeův polynom v Newtonově tvaru.

106

i 0 1 2 3

xi f(xi ), k = 0 −2 −9 0 1 1 3 2 11

f(xi , xi+1 ), k = 1

f(xi , xi+1 , xi+2 ), k = 2

f(xi , . . . , xi+3 ), k = 3

(1 − −9)/(0 − −2) = 5 (3 − 1)/(1 − 0) = 2 (11 − 3)/(2 − 1) = 8

(2 − 5)/(1 − −2) = −1 (8 − 2)/(1 − −2) = 3

(3 − −1)/(2 − −2) = 1

Následuje výpočet Lagrangeova interpolačního polynomu v Newtonově tvaru sage: x sage: sage: x^3 +

var('x') P = -9 + 5*(x+2) - 1*(x+2)*(x - 0) + 1*(x+2)*(x-1)*(x-0) P.expand() x + 1

a ověření správnosti získaného výsledku za pomoci funkce implementované v Sage. sage: R = PolynomialRing(QQ, 'x') sage: f = R.lagrange_polynomial([(-2,-9),(0,1),(1,3),(2,11)]); f x^3 + x + 1

PŘÍKLAD 40 Jsou dány hodnoty napětí naměřené v různých časech. Proložte naměřená data Lagrangeovým polynomem v Newtonově tvaru. Níže uvedený kód obsahuje na svém začátku generování na” měřených“ dat. sage: Uv = 230 sage: pi = 3.14 sage: f = 1 sage: t = [j*0.125 for j in range(20)] sage: u=[Uv*cos(2*pi*f*t[j]) for j in range(20)] sage: o_pair = [(t[j],u[j]) for j in range(20)] sage: R = PolynomialRing(QQ,'t') sage: fun = R.lagrange_polynomial(o_pair) sage: p_lagrange = (plot(fun,0,2, color="red", legend_label='Lagrange', ... title='Lagrangeova interpolace', axes_labels=['t[s]','U[V]']) + ... point(o_pair, legend_label='mereni')) sage: p_lagrange.show() sage: p_lagrange.save("Lagrange.png") sage: print fun 1.54821278562382*t^19 - 36.7908009795200*t^18 + 393.075264439860*t^17 - ... 2487.42652421084*t^16 + 10365.2309043213*t^15 - 30096.7109843768*t^14 + ... 63662.5929060123*t^13 - 103443.634171317*t^12 + 136486.168616036*t^11 - ... 145943.754901973*t^10 + 113457.113109473*t^9 - 58663.9481978659*t^8 + ... 36065.2735467410*t^7 - 33339.1061014539*t^6 + 3903.17895947105*t^5 + ... 14113.4727149469*t^4 + 107.422532226976*t^3 - 4544.00676687586*t^2 + ... 0.300515789687779*t + 230.000000000000

Z grafických výstupů modelu na obrázku 46 je zřejmé, že na výsledky interpolace lze spoléhat pouze v intervalu mezi minimem a maximem z nezávisle proměnné (v příkladu čas t). Mimo tento interval může interpolační polynom nabývat zcela nereálných hodnot neodpovídajících simulovanému jevu.

107

(a) Lagrangeova interpolace

(b) Chování interpolace

Obrázek 46: Ukázka interpolace lagrangeovým polynomem stupně 19.

Splinovací funkce Prokládání naměřených dat polynomem vysokého stupně může vést k získání funkce vykazující nepřirozené oscilace. Tomuto jevu je možné předejít, když setříděná data rozdělíme do menších skupin a interpolaci provádíme po částech polynomy nižšího stupně. Splinovací funkce v Sage je podle všeho kubická.

PŘÍKLAD 41 Uvažujme stejný scénář jako v předchozím příkladu. Jsou dány hodnoty napětí naměřené v různých časech. Proložte naměřená data splinovou (splajnovací) funkcí. Uvedený kód opět začíná vygenerováním naměřených“ dat. ” sage: Uv = 230 sage: pi = 3.14

108

sage: t = [j*0.125 for j in range(20)] sage: f = 1 sage: u=[Uv*cos(2*pi*f*t[j]) for j in range(20)] sage: o_pair = [(t[j],u[j]) for j in range(20)] sage: S = spline(o_pair) sage: p_spline = (plot(S,0,2, color="red", legend_label='spline', ... title='Splinovaci funkce', axes_labels=['t[s]','U[V]']) + ... list_plot(o_pair, plotjoined= False, color='blue')) sage: p_spline.show() sage: p_spline.save("spline.png")

Obrázek 47: Interpolace měření splinovací funkcí. Pro námi uvažovaný případ se výsledky Lagrangeovy interpolace shodují s výstupy využívajícími splajnovací funkce (úsudek na základě grafického srovnání, viz Obrázek 46). To platí pro oblast naměřených dat. Mimo tuto oblast nemá smysl splinovací funkce extrapolovat, protože v takové extrapolaci by byly zohledněny pouze 4 poslední měření (pro kubickou splajnovací funkci).

7.4.2 Regrese Pokud připustíme fakt, že hodnoty získané měřením nejsou stanoveny přesně, pak je vhodnější využít namísto interpolace takzvanou regresi. V případě použití regrese se hledá přibližný, formálně jednoduchý, nikoliv přesný vztah mezi měřenými veličinami. Kritéria pro nalezení takového vztahu mohou vycházet například z minimalizace vybrané (tzv. účelové) funkce. Na základě minimalizace účelové funkce je založena například Metoda nejmenších čtverců (LS ” Least Squares)“. Na jiném principu funguje Úplná metoda nejmenších čtverců (TLS - Total ” 109

Least Squares)“. Využití obou těchto metod je níže ukázáno na příkladu zpracování dat z nepřímého měření elektrického odporu. Metoda nejmenších čtverců Metoda nejmenších čtverců operuje s minimalizací účelové funkce S=

n ∑

(yi − f (xi ))2

i=1

To znamená, že naměřená 2D-data jsou aproximována takovou křivkou, která minimalizuje součet vzdáleností bodů od křivky ve svislém směru (směr osy y). Pokud prokládáme měřenými daty přímku, pak je postup minimalizace následující S(a, b) = ∂S(...) ∂a

=

∂S(...) ∂b

=

n ∑

i=1 n ∑ i=1 n ∑

(yi − axi + b)2 =

n ∑ i=1

yi2 − 2axi yi − 2byi + 2abxi + a2 x2i + b2

2ax2i − 2xi yi + 2bxi = 0 2b − 2yi + 2axi = 0

i=1

Vyřešením soustavy složené z parciálních derivací účelové funkce S(. . .) získáme parametry a, b, které definují hledanou přímku s předpisem y = ax + b. Je to snadné, protože se jedná o řešení dvou lineárních rovnic se dvěma neznámými.

PŘÍKLAD 42 Hodnota nepřímo měřeného elektrického odporu je počítána z naměřených hodnot elektrického napětí a elektrického proudu. Zapojení měřících zařízení ukazuje následující obrázek.

Obrázek 48: Nepřímé měření odporu rezistoru. Při pokusu sestaveného podle uvedeného schématu byly naměřeny hodnoty uvedené v následující tabulce50 . 50

Ve skutečnosti byla data opět generována v Sage.

110

U [V], xi 0,04573733447040 0,37206787799980 0,68273306953338 0,89007675278642 0,96068959920091 1,13171737443341 1,44194258112447 1,50126014907207 1,69500467148787 1,77821243901811 1,80165898386572 2,34442517273023 2,46742838180979 2,74602837003026 3,37940138324154 3,69661673017026 3,85450475629838 3,99274789434560 4,61548761972429 4,73355292742257

Rovnice ∂S(...) ∂a

=

∂S(...) ∂b

=

n ∑ i=1 n ∑

I [mA], yi 0,01374791839983 0,11597655158118 0,20624130974631 0,26972022811710 0,29441665091960 0,34184184377699 0,44562669444350 0,46127700986234 0,52943204559383 0,52070894911111 0,54439192053142 0,69156432228178 0,76159346555334 0,84537476146310 1,05099445651401 1,10188231901274 1,20488121673353 1,22855550608557 1,37513339536447 1,45964061185571

ax2i − xi yi + bxi = 0 b − yi + axi = 0

i=1

mají pro tento případ podobu 135,935656978936 · a + 44,1312940687654 · b + 0,0414861851873037 = 0,0 44,1312940687654 · a + b − 0,0134630011769475 = 0,0 a výsledná regresní přímka má předpis y = 0,000305057773070555 · x + 4,06885608155983 · 10−7 . Odhad velikosti elektrického odporu je převrácenou hodnotou směrnice extrapolované přímky, tedy R = 3278,07 Ω. Grafický výstup modelu bude prezentován níže spolu s výsledkem využití Úplné metody nejmenších čtverců“. ” Metodu nejmenších čtverců je možné aplikovat i na prokládání dat polynomy vyšších stupňů než 1. Nemusíme prokládat data jenom přímkou. Úplná metoda nejmenších čtverců Potřebnou teorii o metodě nejmenších čtverců lze dohledat pod hesly Total Least Squares ” - TLS“, Principal Component Analysis - PCA“ (statistický přístup) nebo Singular Value ” ” Decomposition - SVD“ (algebraický přístup) a zde se jí nebudeme věnovat. Zjednodušeně řečeno, úplná metoda nejmenších čtverců hledá takovou přímku aby součet vzdáleností bodů reprezentujících měření byl od hledané přímky co nejmenší nejen ve směru osy y, ale celkem. Předpokládejme uložení hodnot naměřených napětí a elektrických proudů pořadě v polích u_mer a i_mer. Podobně jako u metody nejmenších čtverců chceme pro nalezení elektrického odporu určit nejprve směrnici extrapolační přímky. Postup jejího nalezení má několik kroků a vysvětluje níže uvedený kód. 1. Od každého prvku pole u_mer je odečtena průměrná hodnota ze všech naměřených napětí. 111

2. Od každého prvku pole i_mer je odečtena průměrná hodnota ze všech naměřených proudů. Odečtení průměrných hodnot vyvolá přesun středu grafické (resp. vektorové) reprezentace měření do počátku souřadné soustavy. 3. Centrovaná“ data jsou uspořádána do matice B. ” 4. Je vypočtena kovarianční matice W B, která přísluší matici posunutých měření B. 5. Jsou nalezena vlastní čísla a vlastní vektory kovarianční matice W B. 6. Ze souřadnic vlastního vektoru příslušejícího v absolutní hodnotě většímu vlastnímu číslu už jde vypočítat směrnici extrapolační přímky a z ní odhad velikosti hledaného odporu. sage: uB = [u_mer[j]-mean(u_mer) for j in range(20)] sage: iB = [i_mer[j]-mean(i_mer) for j in range(20)] sage: B = matrix([(uB[j],iB[j]) for j in range(20)]).transpose() sage: WB = matrix(QQ,B*B.transpose()) sage: print WB.eigenvalues() [4.759039113421180?e-9, 38.55710476829158?] sage: print WB.right_eigenmatrix() ([4.759039113421180?e-9 0] [ 0 38.55710476829158?], [ 1 1] [-3273.320315296367? 0.000305500196644?])

Významnější z vlastních čísel kovarianční matice má hodnotu λ = 38,55710476829158 a příslu→ ší mu vlastní vektor − v = (1; 0,000305500196644)T . Směrnice extrapolační přímky je podílem složek tohoto vektoru a hledaný odpor je její převrácenou hodnotou. Vychází R = 3273, 32031529034 Ω. Ještě snažší cestou k nalezení hodnoty hledaného elektrického odporu je využití funkcí předdefinovaných k prokládání naměřených dat v programu Sage. V kódu je využito načtení uložených dat z předchozích sezení. sage: var('a,b,t') (a, b, t) sage: lsqua = 0.000305057773070555 * t + 4.06885608155983e-7 # data ulozena pomoci save_session(os.path.join("./",'regrese')) sage: load_session(os.path.join("./",'regrese')) sage: model(t) = a*t + b sage: C = matrix(RR,[u_mer, i_mer]).transpose() sage: fit = find_fit(C, model, solution_dict=True) sage: model.subs(fit) t |--> 0.00030550019660604096*t - 9.558918772399863e-07 sage: odporne_ctverce = (plot(model.subs(fit), (t,0,5), color="red", ... legend_label='TLS', title='Voltamperova charakteristika', ... axes_labels=['$U[V]$','$I[A]$']) + point(o_pair_mer, legend_label='mereni') ... + plot(lsqua, (t,0,5), color="green", legend_label = 'LS')) sage: odporne_ctverce.show() sage: odporne_ctverce.save("odporne_ctverce.png")

Výsledek porovnnání je zachycen na Obrázku 49. Z výrazné podobnosti směrnice extrapolované přímky získané v Sage s výsledkem dosaženým vlastním výpočtem (shoda na 6 desetiných míst) lze usuzovat, že v Sage je implementována úplná metoda nejmenších čtverců.

112

Obrázek 49: Shoda výsledku simulací různými metodami. Vzhledem k tomu, že závislost mezi proudem a napětím je lineární tak jsou výsledky využití LS a TLS velmi podobné a regresní přímky se v grafu na Obrázku 49 překrývají. Absolutní členy regresních přímek nevyšly nulové vlivem konečné přesnosti a nastavené tolerance numerických výpočtů. V případě práce s naměřenými daty by se stejným způsobem mohly projevit chyby měření. Výrazněji znatelný rozdíl ve výsledcích obou druhů modelů by byl patrný v případě, kdy bychom regresí aproximovali afinní funkční závislost (tj. přímka by neprocházela počátkem). Modelový příklad afinní závislosti je možné získat prostým přičtením konstanty ke všem hodnotám uloženým v poli i_mer.

7.4.3 Hodnota funkce v bodě Hornerovo schema Nalezení předpisu polynomiální funkční závislosti v naměřených datech má význam pro získávání přibližných hodnot závislé veličiny/proměnné v intevalu měření nezávislé veličiny/proměnné. K určení funkční hodnoty interpolačního polynomu ve vybraném bodě je možné použít nástroj označovaný jako Hornerovo schema“. ” Hornerovo schema je postaveno na rozkladu polynomu P (x) s koeficienty {a0 , . . . , an } do níže zapsané podoby, která rekurzivně zvyšuje stupeň interpolačního polynomu. P (x) = a0 xn + . . . + an P (x) = ((. . . (a0 x + a1 )x + a2 )x + . . . + an−1 )x + an Kromě určování funkční hodnoty polynomu ve zvoleném bodě se Hornerovo schema využívá také k dělení polynomů P (x)/(koefx · x − a).

113

K praktické realizaci výpočtu funkční hodnoty P (x) v bodě x = a se využívá přehledné tabulky. koefx a0 a1 a2 . . . an b(i−1) · a 0 aa0 ab1 . . . abn−1 bi = (ai + b(i−1) · a)/koefx a0 b1 b2 . . . bn = P (a) Když vyjde funkční hodnota P (a) = 0, pak je x = a kořenem rovnice P (x) = 0 a čísla v posledním řádku tabulky představují koeficienty a0 , b1 , . . . , bn−1 podílu polynomů P (x)/(koefx · x − a).

PŘÍKLAD 43 Vydělte polynom x3 − 6x2 + 11x − 6 polynomem (x − 2). 1 1 −6 11 −6 2 0 6 −8 2 1 −4 3 0 Výsledek dělení zadaných polynomů je x2 − 4x + 3. Hodnota polynomu v bodě 2 je 0. Jestliže vyjde s využitím Hornerova schematu funkční hodnota P (a) ̸= 0, potom tato hodnota představuje zbytek po dělení polynomů P (x)/(koefx · x − a).

PŘÍKLAD 44 Vydělte polynom 4x4 − 6x3 + 3x − 5 polynomem (2x − 1). 3 −5 2 4 −6 0 1 0 2 −2 −1 1 2 −2 −1 1 −4 4 Výsledek dělení zadaných polynomů je 2x3 − 2x2 − x + 1 − 2x−1 . Hodnota polynomu v bodě 1/2 je −4.

PŘÍKLAD 45 V kapitole věnované interpolaci jsme proložením navzorkovaného harmonického signálu získali polynom 19. stupně. P (t) = 1,54821278562382 · t19 − 36,7908009795200 · t18 + 393,075264439860 · t17 − 2487,42652421084 · t16 + 10365,2309043213 · t15 − 30096,7109843768 · t14 + 63662,5929060123 · t13 − 103443,634171317 · t12 + 136486,168616036 · t11 − 145943,754901973 · t10 + 113457,113109473 · t9 − 58663,9481978659 · t8 + 36065,2735467410 · t7 − 33339,1061014539 · t6 + 3903,17895947105 · t5 + 14113,4727149469 · t4 + 107,422532226976 · t3 − 4544,00676687586 · t2 + 0,300515789687779 · t + 230,0 Určete hodnotu interpolačního polynomu v bodě t = 1,5 s.

114

Výpočet hodnoty v bodě za použití Hornerova schematu lze snadno algoritmizovat v tabulkovém procesoru nebo jinde. stupeň 1,00 1,50 stupeň 1,00 1,50 stupeň 1,00 1,50

19 1,55 0,00 1,55

18 -36,79 2,32 -34,47

17 393,08 -51,70 341,37 12 11 -103443,63 136486,17 52758,64 -76027,48 -50684,99 60458,68 6 5 -33339,11 3903,18 25289,66 -12074,17 -8049,44 -8170,99

16 15 14 13 12 -2487,43 10365,23 -30096,71 63662,59 -103443,63 512,06 -2963,05 11103,27 -28490,16 52758,64 -1975,37 7402,18 -18993,44 35172,43 -50684,99 10 9 8 7 6 -145943,75 113457,11 -58663,95 36065,27 -33339,11 90688,03 -82883,59 45860,28 -19205,50 25289,66 -55255,73 30573,52 -12803,67 16859,77 -8049,44 4 3 2 1 0 14113,47 107,42 -4544,01 0,30 230,00 -12256,48 2785,49 4339,36 -306,97 -460,00 1856,99 2892,91 -204,64 -306,66 -230,00

Interpolační polynom P (t) má v bodě t = 1,5 s hodnotu −230,0 V. Správnost výsledku lze ověřit na grafu interpolační funkce na obrázku 50.

Obrázek 50: Proložení naměřených dat Lagrangeovým interpolačním polynomem.

Taylorův polynom Taylorův polynom je nástroj k odhadu funkčních hodnot neznámé funkce f (x) na základě známé funkční hodnoty f (a) v bodě x = a a s pomocí známého, diferencovatelného předpisu derivace f ′ (x) funkce f (x). Přesnost odhadu funkční hodnoty f (x) pomocí Taylorova polynomu se zpravidla výrazně zhoršuje se zvyšováním vzdálenosti |x − a| od bodu a se známou funkční hodnotou f (a). f′ (a) f′′ (a) f(n) (a) f(x) = f(a) + (x − a) + (x − a)2 + · · · + (x − a)n + Rn (x) 1! 2! n! 115

kde Rn (x) značí chybu aproximace vzhledem k přesné funkční hodnotě. Rn (x) jde pouze odhadovat, protože f(x) neznáme.

PŘÍKLAD 46 Určete Taylorovy polynomy řádu k = 4 (řád derivace) pro aproximaci hodnot níže uvedených funkcí v blízkosti bodu a = 1. f1 (x) = cos x, f2 (x) = e−x , f3 /x) = ln(x + 1). Ke snadnému výpočtu Taylorových polynomů byly využity příslušné příkazy programu Sage. sage: reset() sage: var('a k x') (a, k, x) sage: f1 = cos(x) sage: f2 = exp(-x) sage: f3 = log(x+1) sage: k = 4; a = 1; sage: t1 = f1.taylor(x,a,k) sage: t2 = f2.taylor(x,a,k) sage: t3 = f3.taylor(x,a,k) sage: p1 = plot(f1, ((-1)*pi,pi), color='blue', legend_label='cos x', ... title='Tyloruv polynom cos x', axes_labels=['$x$','$y$']) + plot(t1,((-1)*pi,pi), ... color='red', legend_label='Taylor') sage: p1.show() sage: p1.save('cosTaylor.png') sage: p2 = plot(f2, ((-1)*pi,pi), color='blue', legend_label='exp(-x)', ... title='Tyloruv polynom exp(-x)', axes_labels=['$x$','$y$']) + ... plot(t2, ((-1)*pi,pi), color='red', legend_label='Taylor') sage: p2.show() sage: p2.save('expTaylor.png') sage: p3 = plot(f3, (0,2*pi), color='blue', legend_label='log(x+1)', ... title='Tyloruv polynom log(x+1)', axes_labels=['$x$','$y$']) + plot(t3, (0,2*pi), ... color='red', legend_label='Taylor') sage: p3.show() sage: p3.save('logTaylor.png')

Aproximace funkčních hodnot uvažovaných funkcí cos x, e−x , log(x + 1) Taylorovými polynomy ukazují po řadě Obrázky 51, 52 a 53. Uvedené grafy prokázaly zhoršování odhadu hodnoty funkce spolu se vzdalováním se od bodu funkční hodnotou kde je funkční hodnota známa.

116

Obrázek 51: Aproximace funkčních hodnot funkce cos x Taylorovým polynomem.

Obrázek 52: Aproximace funkčních hodnot funkce e−x Taylorovým polynomem.

Obrázek 53: Aproximace funkčních hodnot funkce log(x + 1) Taylorovým polynomem.

117

118

8 Závěr Všichni kdo, jsme se na této publikaci podíleli, jsme si přáli, aby se čtenář dostal až k tomuto závěru. Nebylo naším cílem přinést nové myšlenky v oblasti aplikované matematiky. Chtěli jsme jen naznačit, že jsou jiné možné, netradiční, pohledy na problematiku modelování fyzikálních polí. Záměrně jsme použili příklady, které jsou běžně známé, aby tento jiný úhel pohledu vynikl. Stává se, že matematické nástroje používáme při své práci s jistou samozřejmostí a svou pozornost soustředíme jen na to, aby nám naše experimenty vycházely, tak jak si představujeme. Následně se snažíme k fyzikálním výsledkům přizpůsobit matematický aparát. Většinou to děláme ze dvou důvodů. Tím hlavním by mělo být, že chceme mít jistotu, že náš úsudek je správný a že jsme správně pochopili podstatu fyzikálního děje. Velmi často se však stává v dnešní době, která je plná honby za formálními výsledky, že jednoduché fyzikální principy se přizpůsobují nesmyslně složitému matematickému popisu, jen aby se uměle zvýšila vědecká prestiž publikace. Pro tento účel není náš text určen. Pokud v našich stránkách najde kdokoliv inspiraci pro svoji práci a použije z ní i jeden vzoreček, jsme jako autoři spokojeni, že naše práce nebyla zbytečná.

119

120

Použitá literatura a prameny [1] Bálek, R. et al. Povrchové akustické vlny. Academia, 1986. [2] Bosma, W. – Cannon, J. – Playoust, C. The Magma algebra system. I. The user language. J. Symbolic Comput. 1997, 24, 3-4, s. 235–265. ISSN 0747-7171. doi: 10.1006/jsco. 1996.0125. Dostupné z: . Computational algebra and number theory (London, 1993). [3] Haňka, L. Teorie elektomagnetického pole. SNTL/ALFA, 1975. [4] Hindmarsch, A. C. et al. SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers. ACM Transactions on Mathematical Software. September 2005, 31, 3, s. 363–396. [5] Inan, U. S. – Inan, A. Engineering Electromagnetics. Addison-Wesley, 1999. ISBN 0805344233. [6] Kvasnica, J. Matematický aparát fyziky. Academia, 2004. ISBN 8020006036. [7] MATLAB. version 7.10.0 (R2010a). The MathWorks Inc., 2010. [8] Maxwell, J. C. A Dynamical Theory of the Electromagnetic Field. Wipf & Stock Pub, 1996. ISBN 1579100155. [9] Maxwell, J. C. – Physics. Treatise on Electricity and Magnetism, Vol. 1. Dover Publications, 1954. ISBN 0486606368. [10] Maxwell, J. C. – Physics. Treatise on Electricity and Magnetism, Vol. 2. Dover Publications, 1954. ISBN 0486606376. [11] Mayer, D. Teorie elekromagnetického pole 1. a 2. díl. 2004. [12] Merhaut, J. Teoretické základy elektroakustiky. Academia, 1975. [13] Monagan, M. B. et al. Maple 10 Programming Guide. Maplesoft, 2005. [14] Nekvinda, M. – Šrubař, J. – Vild, J. Praktikum z numerické matematiky. 1973. [15] Nye, F. J. Physical Properties of Crystals. Clarendon Press, 1954. [16] Obrazová, H. et al. Příklady z elektrotechniky a elektroniky. Vydavatelství ČVUT, Zikova 4, 166 36 Praha 6, 2000. [17] Rektorys, K. Přehled užité matematiky. 1. a 2. díl. Prometheus, 1995. [18] Sedlák, B. – Štoll, I. Elektřina a magnetizmus. Academia, 2002. ISBN 8020010041. [19] Stein, W. – others. Sage Mathematics Software (Version 6.1.1). The Sage Development Team, 2015. http://www.sagemath.org. [20] The Sage Development Team, 2015. http://sagemath.org/links-components.html. [21] The Sage Development Team, 2015. http://sagemath.org/tour.html. [22] Škvor, Z. Akustika a elektroakustika. Academia, 2001. ISBN 80-200-0461-0. [23] Wolfram Research, Inc. Mathematica (Version 10.1), 2015.

121

Název Autoři Vydavatel Určeno pro čj. RE 93/15 Vyšlo Počet stran Vydání Tiskárna Číslo publikace

Matematický aparát pro modelování fyzikálních polí v příkladech Ing. Lukáš Zedek, Ph.D.; prof. Ing. RNDr. Miloslav Košek, CSc. doc. Ing. Milan Hokr, Ph.D.; Ing. Petr Rálek, Ph.D. Technická univerzita v Liberci odbornou veřejnost, případně pro studenty doktorských studijních programů na FM schváleno rektorátem TU v Liberci dne 9. 10. 2015 listopad 2015 122 1. Vysokoškolský podnik Liberec, s.r.o., Studentská 1402/2, Liberec 55-093-15

Vydání odborné knihy schválila vědecká redakce TUL. Tato publikace neprošla redakční ani jazykovou úpravou. ISBN 978-80-7494-268-6