Symbolické výpočty v 0.2
Jan Schee
Abstract
Vedle numerických výpočtů umožňují dnešní systémy účinně zacházet s algebraickými výrazy, symbolicky řešit algebraické a
diferenciální rovnice. Zvládají také symbolické derivování a integraci.
Cílem tohoto textu je seznámit studenty s řešenímmatematických a fyzikálních problémů pomocí nejrozšířenějšího
nástroje tohoto typu - Wolfram Mathematica (ve švěte OpenSource můžete najít podobný projekt nazvaný Maxima). Stephen
Wolfram je zakladatelem společnosti Wolfram Research, která tento software vyvíji. Aktuální verze (v roce 2015) je 10.02.
Algebraické výrazy
Při spuštění vám Mathematica zpravidla nabídne prázdný notebook do které lze zadávat a vyhodnocovat matematické výrazy.
Chcete-li například znát výsledek výpočtu 1 + 1 tak v notebooku zapíšte tento výraz a výpočet spustíte kombinací kláves
Shift+Enter. V notebooku bude tato operace vypadat takto
In[1]:= 1+1
Out[1]= 2
Mathematica rozlišuje následující matematické operace
Součet x+y+z
Rozdíl x-y
Součin x*y*z nebo xyz
Mocnina x^y
Vyzkoušejte si ve svém notebooku následující příklady:
Sestavování výpočetních operací
Využití předchozích výsledků
V Mathematica se pod výrazem % skrývá předchozí výsledek. Způsoby jak se odkazovat na předchozí výsledky ukazuje tabulka
poslední vygenerovaný výsledek %
předposlední vygenerovaný výsledek %%
k-tý předchozí výsledek %%...% (k-krát)
výsledek na n-tém výstupu Out[n] %n
Vyzkoušejte následující příklady:
In[62]:= 77^2
Out[62]= 5929
In[63]:= %+1
Out[63]= 5930
In[64]:= 3% + %^2+%%
Out[64]= 35188619
In[65]:= %6 + %7
Out[65]= InterpolatingFunction Domain: {{0., 40.}}
Output: scalar
+
θ→ InterpolatingFunction Domain: {{0., 40.}}
Output: scalar
 ,
InterpolatingFunction Domain: {{0., 40.}}
Output: scalar
+
θ′ → InterpolatingFunction Domain: {{0., 40.}}
Output: scalar
 
Definice proměnných
Když provádíte dlouhé matematické operace, pak užitečné jednotlivé mezivýsledky pojmenovat. Tak jako v matematice nebo
programovacích jazycích i zde k tomu účelu slouží koncept proměnné.
Následjící příklady nejlépe výsvětlí jak s proměnnými v Mathematica zacházet.
◼ Nastavení proměnné x na hodnotu 5
In[66]:= x=5
Out[66]= 5
◼ Kdykoliv se teď objeví ve výrazech x, Mathematica na jeho místo dosadí číslo5
In[67]:= x^2
Out[67]= 25
◼ Chceme-li do proměnné x vložit jinou hodnotu, tak samozřejmě můžeme
In[68]:= x=4+7
Out[68]= 11
◼ Nastavme proměnnou pi na hodnotu číslaπ s přesností na 40 desetinných míst
In[69]:= pi=N[π, 40]
Out[69]= 3.141592653589793238462643383279502884197
In[70]:= pi
Out[70]= 3.141592653589793238462643383279502884197
◼ Druhá mocnina číslauloženého v proměnné pi pak bude
In[71]:= pi^2
Out[71]= 9.86960440108935861883449099987615113531
Následující tabulka shrnuje způsoby definice proměnných.
přiřazení hodnoty promměné x x = hodnota
přiřazení hodnoty proměnným x a y x=y=hodnota
odstranění jakékoliv hodnoty přiřazené x x=. nebo Clear[x]
2 symbolicke_vypocty.nb
Seznamy objektů
Během býpočtů je občas výhodné seskupit dohromady několik objektů a pracovat s nimi jako s jednou entitou. Systém
Mathematica pro tento účel nabízí koncept seznamu - List.
Například seznam {2,3,5} je kolekcí tří objektů se kterým můžeme, v řadě algebraických manipulací, zacházet jako s jedním
objektem. Vyzkoušejte následující příklady.
◼ Seznam tří čísel
In[72]:= {2, 3, 5}
Out[72]= {2, 3, 5}
◼ Umocnění každého číslav seznamu a přidání jedničky
In[73]:= {2, 3, 5}^2+1
Out[73]= {5, 10, 26}
◼ Rozdíl dvou seznamů
In[74]:= {9, 4, 3}-{3, 2, 5}
Out[74]= {6, 2, -2}
◼ Hodnota předchozího seznamu
In[75]:= %
Out[75]= {6, 2, -2}
◼ Aplikace matematických funkcí na seznam
In[76]:= Exp[%]//N
Out[76]= {403.429, 7.38906, 0.135335}
◼ Přiřazení proměnné v seznam
In[77]:= v={2, 4, 5}
Out[77]= {2, 4, 5}
◼ Kdekoliv se teď v matematických operacích v objeví tak bude nahrazenou seznamem
In[78]:= v/(v-1)
Out[78]= 2,
4
3
,
5
4

Manipulace s prvky seznamu
Na prvky seznamu se odkazujeme prostřednictvím jeho indexu. Číslováníindexu začíná od 1. V následující tabulce jsou shrnuty
základní informace o manipulaci se seznamy.
seznam {a,b,c}
i-tý prvek seznamu Part[seznam, i] nebo seznam[[i]]
seznam i-tých a j-tých prvků seznamu Part[seznam,{i,j}] nebo seznam[[{i,j}]]
Vyzkoušejte následující příklady.
◼ Získání třetího prvku seznamu
In[79]:= {2, 9, 7, 5, 4}[[3]]
Out[79]= 7
◼ Získání seznamu prvků seznamu
In[80]:= {2, 9, 7, 5, 4}[[{4, 1, 5}]]
Out[80]= {5, 2, 4}
symbolicke_vypocty.nb 3
◼ Přiřazení seznamu proměnné v
In[81]:= v={7, 9, 4}
Out[81]= {7, 9, 4}
◼ Získání druhého prvku zeznamu proměnné v
In[82]:= v[[2]]
Out[82]= 9
Čtyři typy závorek v Mathematica
V předchozím jsme se setkali se třemi typi závorek. Každá z nich má specifický význam a nelze je zaměňovat. Následující
tabulka shrnuje tyto čtyřitypy závorek.
kulaté závorky pro grupování (výraz)
hranaté závorky pro funkce f[x]
kroucené závorky pro seznamy {a,b,c}
dvojité závorky pro indexaci v[[i]]
Posloupnost matematických operací
Matematické operace zpravidla provádíme v několika krocích. Pokud chceme můžeme každý krok zapsat na zvláštní řádek.
Občas se ale ukáže užitečné seskupit jednotlivé kroky do jednoho řádku. V takovém případě jsou jednotlivé kroky od sebe
odděleny středníkem.
proveď několik operací a zobraz výsledek poslední expr1;expr2;expr3
proveď operace ale výsledky nedávej na výstup expr1;expr2;
◼ Budou provedeny tři operace a výsledek poslední je zobrazen ve výstupu
In[83]:= x=4;y=6;z=y+6
Out[83]= 12
◼ Vložením středníku na konec řádku říkámesystému Mathematica aby výstup nezobrazoval
In[84]:= x=67-5;
◼ Slále můžeme použít % k zobrazení výsledku
In[85]:= %
Out[85]= 62
Moduly a bloky
Mathematica umožňuje určitou sekvenci příkazů sloučit do jednoho prostředí ve které lze používat lokálně definované
proměnné. K tomuto slouží Module. Vyzkoušejme následující příklad
◼ Definice modulu s lokálními proměnnými u,v a asociace s funkcí f
In[86]:= f[x_]:=Module[{u, v},
u=x^2+1;
v=u^2;
v
];
◼ Nastavení globálních proměnných u,v
In[87]:= u=1;v=2;
◼ Zavolání funkce f a zobrazení hodnot globálních proměnných u,v
4 symbolicke_vypocty.nb
In[88]:= f[5]
u
v
Out[88]= 676
Out[89]= 1
Out[90]= 2
Další nástroj umožnující sloučení sekvence příkazů do jednoho bloku je Block. Tato struktura funguje tak, že globální
proměnné se při spuštění Blocku nejprve vyprázdní, v rámci Blocku se naplní novými čísly a po ukončení Blocku jsou globální
proměnné opět nastaveny na původní hodnoty.
◼ Nejprve se nastavi hodnoty proměnných a,b,c. V rámci struktury Block se změní jejich hodnoty. Na konci této struktury
jsou obnoveny původní hodnoty a,b,c
In[91]:= a=1;b=2;c=3;
In[92]:= Block[{a, b, c},
a=4;
b=5;
c=a+b
]
Out[92]= 9
In[93]:= a
b
c
Out[93]= 1
Out[94]= 2
Out[95]= 3
In[96]:=
Definice funkce v systému Mathematica
Funkce je v systému Mathematica jednoznačně identifikována hranatými závorkami, [].
◼ Definujme funkci F s jedním argumentem x jako x^2
In[2]:= F[x_]=x^2
Out[2]= x2
◼ Zpožděné přiřazení (Výraz na pravé straně je vyhodnocen až při použití funkce G)
In[3]:= G[x_]:=x^2
Algebraické výpočty - symbolické výpočty
Rozdíl mezi numerickým a symbolickým výpočtem ukazují následující ukázky.
◼ Numerický výpočet
In[99]:= 13+2-5
Out[99]= 10
◼ Symbolický výpočet
symbolicke_vypocty.nb 5
In[100]:= 2x+1
Out[100]= 125
K symbolické manipulaci s algebraickými výrazy Mathematica definuje celou řadu funkcí. Některé z nich zde představíme v
rámci příkladů. Zbylé funkce snadno naleznete v nápovědě systému Mathematica.
◼ Expand[] - tato funkce roznásobuje a umocňuje algebraické výrazy
In[101]:= (x+2y+1) (x-2)^2
Out[101]= 270000
In[102]:= Expand[%]
Out[102]= 270000
◼ Factor[] - tato funkce dělá pravý opak funkce Expand[]
In[103]:= Factor[%]
Out[103]= 270000
In[104]:=
Mathematica také nabízí funkce, které se snaží standartními algebraickými operacemi výrazy zjednodušit.
◼ Simplify[] - tato funkce se snaží zjednodušit algrebraický výraz použitím standartních algebraických transformací
In[4]:= Simplify [x^2-2x+1]
Out[4]= (-1+x)2
In[5]:= Integrate[1/(x^4-1), x]
Out[5]= -
ArcTan[x]
2
+
1
4
Log[1-x]-
1
4
Log[1+x]
In[6]:= D[%, x]
Out[6]= -
1
4(1-x)
-
1
4(1+x)
-
1
2(1+x2)
In[7]:= Simplify [%]
Out[7]=
1
-1+x4
◼ FullSimplify[] - tato funkce se snaží zjednodušit algebraický výraz zapojením širšího spektra transformací
In[8]:= Simplify [Gamma [x] Gamma [1-x]]
Out[8]= Gamma [1-x] Gamma [x]
In[9]:= FullSimplify [Gamma [x] Gamma [1-x]]
Out[9]= π Csc[π x]
Pro algebraickou manipulaci s racionálními Mathematica definuje další sérii užitečných funkcí.
◼ ExpanAll[ ] - tato funkce aplikuje Expand[] všude (třeba , v případě racionálních výrazů, v čitatelii jmenovateli )
In[10]:= e=(x-1)^2(x+2)/((1+x) (x-3)^2)
Out[10]=
(-1+x)2
(2+x)
(-3+x)2
(1+x)
In[11]:= ExpandAll[e]
Out[11]=
2
9+3x-5x2 +x3
-
3x
9+3x-5x2 +x3
+
x3
9+3x-5x2 +x3
6 symbolicke_vypocty.nb
◼ Together[] - tato funkce dá všechny členyna společný jmenovatel
In[12]:= Together[%]
Out[12]=
2-3x+x3
(-3+x)2
(1+x)
◼ Apart[] - v tomto případě rozdělí složidý zlomek na dílčí, jednodušší zlomky
In[13]:= Apart[%]
Out[13]= 1+
5
(-3+x)2
+
19
4(-3+x)
+
1
4(1+x)
◼ Cancel[] - vykrácení společné faktory mezi čitatelma jmenovatelem
In[115]:=
Pokud máme výraz obsahující více proměných pak ,můžeme ve výrazu uspořádat členy do skupin tak, že jedna nebo druhá
proměnná dominuje.
seskupení mocnin proměnné x ve výrazu expr Collect[expr,x]
vytknout členyvýrazu expr, které nezávisí na x FactorTerms[expr,x]
◼ Definujme výraz dvou proměnných
In[14]:= v=Expand[(3+2x)^2(x+2y)^2]
Out[14]= 9x2 +12x3 +4x4 +36xy+48x2 y+16x3 y+36y2 +48xy2 +16x2 y2
◼ Seskupíme členyvýrazu v obsahující stejnou mocninu proměnné x
In[15]:= Collect[v, x]
Out[15]= 4x4 +36y2 +x3 (12+16y)+x2 (9+48y+16y2)+x(36y+48y2)
◼ Seskupíme členyvýrazu v obsahující stejnou mocninu proměnné y
In[16]:= Collect[v, y]
Out[16]= 9x2 +12x3 +4x4 +(36x+48x2 +16x3) y+(36+48x+16x2) y2
◼ Vytknutí členůve výrazu v, které neobsahují proměnnou x
In[17]:= FactorTerms [v, y]
Out[17]= (9+12x+4x2) (x2 +4xy+4y2)
Zatím jsme se zabývali pouze algebraickými výrazy, Mathematica nabízí taky funkce umožňující manipulaci s trigonometrickými
výrazy.
Rozložení trigonometrického výrazu expr na sumu členů TrigExpand[expr]
Vytvoření foktorovaného tvaru trigonometrického výrazu expr TrigFactor[expr]
Redukce trigonometrického výrazu pomocí násobků úhlů TrigReduce[expr]
Převedení trigonometrických funkcí na exponenciální TrigToExp[expr]
Převedení exponenciálních funkcí na trigonometrické ExpToTrig[expr]
◼ Rozložení trigonometrického výrazu
In[18]:= TrigExpand[Tan[x] Cos[2x]]
Out[18]=
3
2
Cos[x] Sin[x]-
Tan[x]
2
-
1
2
Sin[x]2
Tan[x]
◼ Získání faktorovaného tvaru trigonometrického výrazu
symbolicke_vypocty.nb 7
In[19]:= TrigFactor[%]
Out[19]= 2Sin
π
4
-x Sin
π
4
+x Tan[x]
◼ Redukce výrazu pomocí násobků úhlů
In[20]:= TrigReduce[%]
Out[20]= -
1
2
Sec[x] (Sin[x]-Sin[3x])
Zjednodušení výrazů pomocí “Assumptions”
Využijme funkci Simplify ve tvaru Simplify[expr,assum], kde výraz expr je zjednodušen na základě předpokladů assum.
◼ Mathematica automaticky nezjednošuje obecný výraz x2
protože tento výraz neplatí pro všechna x
In[123]:= Simplify [Sqrt[x^2]]
Out[123]= 62
◼ Řekneme Mathematica aby předpokládala, že platí x>0
In[21]:= Simplify [Sqrt[x^2], x>0]
Out[21]= x
◼ Další výraz Mathematica taky automaticky nezjednoduší
In[22]:= 2a+2Sqrt[a-Sqrt[-b]] Sqrt[a+Sqrt[-b]]
Out[22]= 2a+2 a- -b a+ -b
◼ Za předpokladu, že je a>0 a současně b>0 pak Mathematica výraz zjednoduší
In[23]:= Simplify [%, a>0&&b>0]
Out[23]= 2a+ a2 +b 
V dalším využijeme funkci Element[x,dom], jejíž výstupem je tvrzení, že x je prvkem množiny (intervalu) dom.
◼ Necháme opět Mathematica upravit výraz x2
za předpokladu, že x je reálné číslo
In[24]:= Simplify [Sqrt[x^2], Element [x, Reals]]
Out[24]= Abs[x]
◼ V další ukázce zjednodušíme goniometrickou funkci, za předpokladu, že proměnná i je celočíselná
In[25]:= Simplify [Sin[x+2nPi], Element [n, Integers]]
Out[25]= Sin[x]
Části algebraických výrazů
Mathematica vám umožní prostřednictvím integrovaných funkcí vybírat z výrazu jeho části, které lze dál použít k následujícím
výpočtům. Zde prezentujeme pět těcho funkcí.
Coefficient [expr, form ] vrátí z výrazu expr koeficient v požadovaném tvaru form
Exponent[expt, form] vrátí z výrazu expr maximální mocninu ve tvaru form
Part[expr,n] nebo expr[[n]] vrátí n-tý členve výrazu expr
Numerator[expr] vrátí čitatelvýrazu expr
Denominator[expr] vrátí jmenovatel výrazu expr
Na následujících příkladech si vyzkoušejte požití těchto funkcí.
8 symbolicke_vypocty.nb
◼ Uvažujme algebraický výraz
In[26]:= e=Expand[(1+3x+4y^2)^2]
Out[26]= 1+6x+9x2 +8y2 +24xy2 +16y4
◼ Koeficient u x ve výrazu e bude
In[27]:= Coefficient[e, x]
Out[27]= 6+24y2
◼ Nejvyšší mocnina y ve výrazu e
In[28]:= Exponent[e, y]
Out[28]= 4
◼ Čtvrtý členve výrazu e
In[29]:= Part[e, 4]
Out[29]= 8y2
In[30]:= e[[4]]
Out[30]= 8y2
◼ Nyní uvažujme racionální výraz
In[31]:= r=(1+x)/(2(2-y))
Out[31]=
1+x
2(2-y)
◼ Čitatel výrazu r bude
In[32]:= Numerator [r]
Out[32]= 1+x
◼ Jmenovatel výrazu r bude
In[33]:= Denominator [r]
Out[33]= 2(2-y)
Základní matematické operace pro symbolické výpočty
Následující tabulka obsahuje seznam funkcí, které Mathematica poskytuje pro diferenciální počet.
parciální derivace ∂f
∂x
D[f,x]
neurčitý integrál ∫f ⅆ x Integrate[f,x]
suma ∑i =imin
imax
f Sum[f,{i,imin,imax}]
řešení rovnice pro proměnnou x Solve[lhs==rhs,x]
rozvoj funkce f do mocninné řadyv okolí x=x0 Series[f,{x,x0, order}]
limita lim x->x0f(x) Limit[x,x->x0]
Dataminimalizace funkce vzhledem k x Minimize[f,x]
◼ derivace funkce xn
In[34]:= D[x^n, x]
Out[34]= nx-1+n
◼ třetí derivace funkce xn
symbolicke_vypocty.nb 9
In[35]:= D[x^n, {x, 3}]
Out[35]= (-2+n) (-1+n) nx-3+n
Mathematica kromě parciálních derivací nabízí i totální derivace. K tomu definuje funkci Dt[f,x].
◼ Totální derivace funkce xn
. Mathematica implicitně předpokládá, že n je funkcí x.
In[36]:= Dt[x^n, x]
Out[36]= xn 
n
x
+Dt[n, x] Log[x]
Zde výraz Dt[n,x] znamená dn
dx
.
◼ Integrací funkce xn
obdržíme
In[37]:= Integrate[x^n, x]
Out[37]=
x1+n
1+n
◼ Nebo integrací složitější funkce
In[38]:= Integrate[1/(x^4-a^4), x]
Out[38]= -
ArcTanx
a

2a3
+
Log[a-x]
4a3
-
Log[a+x]
4a3
◼ Určitý integrál je určený mezemi. Zkuste si následující příklad ∫a
b
sin 2
(x) ⅆ x
In[39]:= Integrate[Sin[x]^2, {x, a, b}]
Out[39]=
1
2
(-a+b+Cos[a] Sin[a]-Cos[b] Sin[b])
◼ Mathematica nám nedá symbolický výraz pro určitý integrál ∫0
1
xx
ⅆ x
In[40]:= Integrate[x^x, {x, 0, 1}]
Out[40]= 
0
1
xx ⅆ x
◼ Nicméně pořád můžete dostat numerický výsledek
In[41]:= N[%]
Out[41]= 0.783431
In[145]:=
Rovnice
Operace “=” je operace přiřazení. Pro rovnost je v Mathematica vyhrazena operace “==”. Vyzkoušejme následující příklady.
◼ Přiřadíme proměnné a číslo4
In[42]:= a=4
Out[42]= 4
◼ Když se zeptáte co je v a získáte 4, tj.
In[43]:= a
Out[43]= 4
◼ Můžeme taky otestovat jaká hodnota je v a
In[44]:= a⩵4
Out[44]= True
10 symbolicke_vypocty.nb
In[45]:= a⩵7
Out[45]= False
◼ Vyčistíme proměnnou a
In[46]:= a=.
◼ Když teď budeme chtít otestovat je -li v a třeba číslo5 pak Mathematica se nebude ani snažit najít odpověď, protože a
nemá definoanou hodnotu.
In[47]:= a⩵5
Out[47]= a⩵ 5
◼ A když nahradíme a specifickou hodnotou, třeba 4 dostaneme výsledek předchozí akce false.
In[48]:= %/.a→4
Out[48]= False
◼ Pokud budou dva symbolické výrazy identické pak Mahematica vrátí true
In[49]:= 2x+x^2⩵2x+x^2
Out[49]= True
◼ Ovšem stačí upravit výraz na pravé straně (například) a Mathematica se nebude ani snažit najít odpověď.
In[50]:= 2x+x^2⩵x(2+x)
Out[50]= 2x+x2 ⩵ x(2+x)
◼ Mějme nyní rovnici
In[51]:= x^2-3x+5⩵0
Out[51]= 5-3x+x2 ⩵ 0
◼ Přiřadíme jí jméno
In[52]:= eq=%
Out[52]= 5-3x+x2 ⩵ 0
In[53]:= eq
Out[53]= 5-3x+x2 ⩵ 0
Pro řešeníalgebraických rovnic, Mathematica poskytuje funkci Solve.
Solve[lhs == rhs, x]
Řeší rovnici vzhedem k x, dostáváme toto řešeníve formě
pravidel vzhledem k x
x/.solution ze seznamu řešenízíkáme jednotlivá řešení
expr/.solution
použití seznamu řešení(pravidel) k zístání hodnoty pro výraz
expr
◼ Mathematica umí vždy vyřešit agebraickou rovnici jedné proměnné symbolicky, pokud je její řádmenší než 5
symbolicke_vypocty.nb 11
In[54]:= Solve[x^4-2x^3+3x^2-6x+12⩵0, x]
Out[54]= x→
1
2
-
1
2
-1+
13
-11+2ⅈ 519 
1/3
+-11+2ⅈ 519 
1/3
-
1
2
-2-
13
-11+2ⅈ 519 
1/3
--11+2ⅈ 519 
1/3
-
8
-1+ 13
-11+2ⅈ 519 
13
+-11+2ⅈ 519 
1/3
,
x→
1
2
-
1
2
-1+
13
-11+2ⅈ 519 
1/3
+-11+2ⅈ 519 
1/3
+
1
2
-2-
13
-11+2ⅈ 519 
1/3
--11+2ⅈ 519 
1/3
-
8
-1+ 13
-11+2ⅈ 519 
13
+-11+2ⅈ 519 
1/3
,
x→
1
2
+
1
2
-1+
13
-11+2ⅈ 519 
1/3
+-11+2ⅈ 519 
1/3
-
1
2
-2-
13
-11+2ⅈ 519 
1/3
--11+2ⅈ 519 
1/3
+
8
-1+ 13
-11+2ⅈ 519 
13
+-11+2ⅈ 519 
1/3
,
x→
1
2
+
1
2
-1+
13
-11+2ⅈ 519 
1/3
+-11+2ⅈ 519 
1/3
+
1
2
-2-
13
-11+2ⅈ 519 
1/3
--11+2ⅈ 519 
1/3
+
8
-1+ 13
-11+2ⅈ 519 
13
+-11+2ⅈ 519 
1/3

◼ První kořen předchozí rovnice pak bude
In[55]:= x/.%[[1]]
Out[55]=
1
2
-
1
2
-1+
13
-11+2ⅈ 519 
1/3
+-11+2ⅈ 519 
1/3
-
1
2
-2-
13
-11+2ⅈ 519 
1/3
--11+2ⅈ 519 
1/3
-
8
-1+ 13
-11+2ⅈ 519 
13
+-11+2ⅈ 519 
1/3
12 symbolicke_vypocty.nb
In[56]:= N[%]
Out[56]= -0.611432-1.69854ⅈ
◼ Dokáže si poradit i s některými rovnicemi vyššího řádunež 4
In[57]:= Solve[x^6⩵1, x]
Out[57]= {x→ -1}, {x→ 1}, x→-(-1)1/3
, x→ (-1)1/3
, x→ -(-1)2/3
, x→ (-1)2/3

◼ Ovšem existují algebraické rovnice jejichž řešeníneumí Mathematica vyjádřit symbolicky
In[58]:= Solve[x^5-3x^4+6x^3-2x+16⩵0, x]
Out[58]= x→Root16-2#1+6#13
-3#14
+#15
&, 1,
x→ Root16-2#1+6#13
-3#14
+#15
&, 2, x→Root16-2#1+6#13
-3#14
+#15
&, 3,
x→ Root16-2#1+6#13
-3#14
+#15
&, 4, x→Root16-2#1+6#13
-3#14
+#15
&, 5
◼ V takovém případě je nutné sáhnout k numerickému řešení
In[59]:= N[%]
Out[59]= {{x→-1.18751}, {x→0.458263-1.48994ⅈ},
{x→0.458263+1.48994ⅈ}, {x→1.63549-1.69411ⅈ}, {x→ 1.63549+1.69411ⅈ}}
◼ V některých případech dokáže Mathematica řešiti rovnice zahrnující i jiné funkce
In[60]:= Solve[Sin[x]⩵a, x]
Out[60]= {{x→ConditionalExpression[π-ArcSin[a]+2π C[1], C[1]∈Integers]},
{x→ConditionalExpression[ArcSin[a]+2π C[1], C[1]∈Integers]}}
◼ Ve většině nealgebraických rovnic ovšem neuspěje
In[61]:= Solve[Cos[x]⩵x, x]
Out[61]= Solve[Cos[x]⩵ x, x]
◼ Pak je možné hledat přibližné řešenínumericky
In[62]:= FindRoot[Cos[x]⩵x, {x, 0}]
Out[62]= {x→0.739085}
Mathematica samozřejmě umí rešit systémy rovnic pro více neznámých. V takovém případě se používá ve formátu jak uvádí
tabulka
Solve[{lhs1 == rhs1, lhs2 == rhs2, ...}, {z, y, ...}] řeší systém rovnic pro proměnné x,y,...
◼ Řešme například systém dvou lineárních rovnic pro dvě neznámé x,y
In[63]:= Solve[{ax+y⩵0, 2x+(1-a) y⩵1}, {x, y}]
Out[63]= x→ -
1
-2+a-a2
, y→ -
a
2-a+a2

◼ Mathematica dokáže řešitsystém nelineárních algebraických rovnic
In[64]:= Solve[{x^2+y^2⩵1, x+3y⩵0}, {x, y}]
Out[64]= x→ -
3
10
, y→
1
10
, x→
3
10
, y→ -
1
10

◼ Toto řešeníteď můžeme vzít a použít k vyčíslení výrazu (řešení je totiž vráceno ve formě pravidelm , jako x-> 4)
In[65]:= x+y/.%
Out[65]= -
2
5
,
2
5

◼ Když pracujeme se systémem rovnic pak bývá užitečné reorganizovat tento systém eliminací jedné z proměnných
symbolicke_vypocty.nb 13
In[66]:= Eliminate [{ax+y⩵0, 2x+(1-a) y⩵1}, y]
Out[66]= (2-a+a2) x⩵ 1
Numerické výpočty v Mathematica
Mathematica implementuje komplexní numerické metody do funkcí umožňující řešit daný problém voláním vhodných funkcí
a předávání parametrů získaných v předchozích výpočtech. Nakonec, zde máme celou řadu funkcí, které nám dovolí získané
výsledky vizualizovat. My se zde omezíme na funkce umožňující numerickou integraci, řešení diferencialních rovnic, řešení
algebraických rovnic, hledání kořene transcendentních rovnic.
NIntegrate[f , {x, xmin , xmax }]
Numericky integruje funkci f vzhledem k x v mezich od xmin
do xmax
NDSolve[{odeys,initial conditions},{x,y},{t,tmin,tmax}]
Numericky řešísystém diferenciálních rovnic pro funkce x,y
vzhedem k parametru t od tmin do tmax
NSolve[lhs==rhs,x] Numericky řešíalgebraickou rovnici pro neznámou x
FindRoot[lhs==rhs,{x,x0}] Hledá kořen rovnici x s iniciálním odhadem x0
In[171]:=
In[172]:=
Aplikace
V této částijsou prezentovány problémy řešenév prostředí Mathematica .
Jednoduché kyvadlo
Budeme hledat numerické řešenírovnice popisující jednoduché kyvadlo. K popisu systému použijeme polární souřadnice (r,θ).
Nechť má kyvadlo délku l a nechť se pohybuje v tíhovém poli g.
14 symbolicke_vypocty.nb
In[173]:=
Out[173]=
Lagrangián této soustavy potom je
L = T - V =
1
2
l θ


2
+ g l Cos(θ). (0.1)
Pohybové rovnice pak určíme z Eulerových-Lagrangeových rovnic, které mají v našem případě tvar
∂L
∂θ
-
d
dt
∂L
∂θ
 = 0 ⇒ θ
..
+
g
l
Sin (θ) = 0. (0.2)
Kyvadlo se pohybuje v konzervativním poli a tak se jeho celková energie zachovává. Velikost této energie udává následující výraz
E = T + V =
1
2
l θ


2
- g l Cos(θ). (0.3)
symbolicke_vypocty.nb 15
Pro kyvadlo uvedené do pohybu s počáteční výchylkou θ0a počáteční rychlostí θ

0=0 bude tato energie
E = -g l Cos[θ0]. (0.4)
◼ Vyčistíme a nastavíme parametry l a g
In[67]:= Clear[l,g]
g=1;
l=5;
◼ Nastavíme počáteční podmínky
In[70]:= θ0=-20.Degree;
◼ Řešíme pohybovou rovnici v časechod t=0 do t=tmax
In[71]:= tmax =40;
dsol=NDSolve[{θ''[t]+ g/l Sin[θ[t]]⩵0,θ[0]⩵θ0,θ'[0]⩵0},{θ,θ'},{t,0,tmax }];
◼ Získáme interpolované funkce řešení
In[73]:= fθ=θ/.dsol[[1]];
fv=θ'/.dsol[[1]];
◼ Zobrazíme výsledek v grafu
In[75]:= Plot[{fθ[t],0},{t,0,tmax }]
Out[75]=
10 20 30 40
-0.3
-0.2
-0.1
0.1
0.2
0.3
◼ Ze znalosti celkové energie určeme periodu kmitu
In[76]:= E0=-g l Cos[θ0]
T=-2 NIntegrate[l/Sqrt[Abs[2(E0+g l Cos[θ])]],{θ,-θ0,θ0}]
Out[76]= -4.69846
Out[77]= 14.1574
protože zřejmě platí
E =
1
2
l θ


2
- g l Cos(θ) ⇒ θ

=
2 (E + g l Cos(θ))
l
⇒ T = 2 
-θ0
θ0 l
2 (E + g l Cos(θ))
ⅆ θ. (0.5)
◼ Ověříme do jaké míry se zachovává energie při integraci
16 symbolicke_vypocty.nb
In[78]:= Plot[(1/2(l fv[t])^2- g l Cos[fθ[t]])-E0,{t,0,tmax }]
Out[78]=
10 20 30 40
-1. ×10-6
-8. ×10-7
-6. ×10-7
-4. ×10-7
-2. ×10-7
2. ×10-7
Během integrace dochází k numerické disipaci energie!
Pohyb testovací částice v centálním, konzervativním poli
Máme testovací částici, pohybující se v gravitačním gravitačním poli tří pevných těles. Celkové pole je popsáni gravitačním
potenciálem
ϕtot = ϕI + ϕII + ϕIII (0.6)
kde je
ϕi = -
GMi
r - r
i
. (0.7)
Protože je pole konzervativní, tak se zachovává celková, specisfická energie E testovací částice
E =
1
2
v2
+ ϕtot. (0.8)
Pohyb budeme sledovat v kartézké soustavě, kde polohu každého objektu popisuje vektor
r = x e
x + y e
y + z e
z (0.9)
a rychost určuje zřejmě vektor
v =
d r
dt
=
dx
dt
e
x +
dy
dt
e
y +
dz
dt
e
z = vx e
x + vy e
y + vz e
z. (0.10)
Pohybová rovnice testovací částicepak zřejmě bude
d v
dt
= -∇ϕtot. (0.11)
◼ Vyčistíme proměnné G (grav. konstanta),M1,M2,M3 (hmotnosti objektů generujících grav. pole),r1,r2,r3 (polohové vektory
těchto těles)
In[105]:= Clear[G,M1,M2,M3,r1,r2,r3];
◼ Vektor celkové gravitační síly v mítě určeném polohovým vektorem r
In[106]:= Ftot[r_?VectorQ]:=-G M1 (r-r1)/((r-r1).(r-r1))^(3/2)-G M2 (r-r2)/((r-r2).(r-r2))^(3/2)-G
◼ Nastavení konfigurace systému tří hmotných center
In[107]:= G=1;M1=0.5;M2=1.5;M3=2.5;
r1={0,5,0};
r2={0,-5,0};
r3={-5,0,0};
◼ Počáteční podmínky testovacé částice
symbolicke_vypocty.nb 17
In[111]:= v0={0,0.1,0};(*pocatecni rychlost*)
r0={20,0,0};(*pocatecni polohovy vektor*)
◼ Integrace pohybove rovnice testovaci castice v case of t=0 do t=200
tmax =200;
dsol=NDSolve[{r''[t]⩵Ftot[r[t]],r'[0]⩵v0,r[0]⩵r0},{r,r'},{t,0,tmax }]
Out[117]= r→ InterpolatingFunction Domain: {{0., 200.}}
Output dimensions: {3}
,
r′ → InterpolatingFunction Domain: {{0., 200.}}
Output dimensions: {3}

◼ Získání interpolovaného řešenípro radiální průvodič testovací částicea vektor rychlosti
In[118]:= fr=r/.dsol[[1]]
fv=r'/.dsol[[1]]
Out[118]= InterpolatingFunction Domain: {{0., 200.}}
Output dimensions: {3}

Out[119]= InterpolatingFunction Domain: {{0., 200.}}
Output dimensions: {3}

◼ Zobrazení trajektorie společně s polohami gravitujících center
In[120]:= pl=ParametricPlot3D [fr[t],{t,0,tmax }];
sp=Graphics3D[{Red,Sphere[r1,1/2],Blue,Sphere[r2,1/2],Green,Sphere[r3,1/2]}];
Show[sp,pl]
Out[122]=
◼ Zachovávající se energie testovací částiceje
In[123]:= Energie0=1/2 v0.v0 -G M1/Sqrt[(r0-r1).(r0-r1)]-G M2/Sqrt[(r0-r2).(r0-r2)]-G M3/Sqrt[(r0-r3
Out[123]= -0.192014
◼ Energie jako funkce
In[124]:= Energie[v_?VectorQ,r_?VectorQ]:=1/2 v.v -G M1/Sqrt[(r-r1).(r-r1)]-G M2/Sqrt[(r-r2).(r-r2)]
◼ Graf funkce Energie[t]-Energie0
18 symbolicke_vypocty.nb
In[125]:= Plot[Energie[fv[t],fr[t]]-Energie0,{t,0,tmax }]
Out[125]=
50 100 150 200
-1.5 ×10-6
-1. ×10-6
-5. ×10-7
5. ×10-7
Rovnice vedení tepla
Uvažujme tenkou tyč délky L a s rozložením teploty T(x,t) v místě x a časet. Udržujme konce tyče na konstantí teplotě
T(0, t) = T(L, t) = T0. (0.12)
Iniciální distribuce teploty je T(x,0), distribuce teploty v tyči v libovolném časeurčuje rovnice
∂T
∂t
= α
∂2
T
∂x2
(0.13)
Dále určeme následující integrál
d
dt

0
1
T(x, t) ⅆ x = 
0
L d
dt
T(x, t) ⅆ x = α 
0
L ∂
∂x
∂T
∂x
ⅆ x = α
∂T
∂x L
- α
∂T
∂x 0
. (0.14)
◼ Nechť naši tyč charakterizuje parameter α a délka tyče L
In[126]:= α=1;
L=10;
◼ Počáteční rozložení teploty v tyči určuje funkce T0
In[128]:= T0[x_]=Sin[x/L Pi];
◼ Hledání numerického řešení,při splnění okrajových podmínek T(0,t)=T(L,t)=0.
In[131]:= T=.
tmax =40;
dsol=NDSolve[{D[T[x,t],t]-α D[T[x,t],{x,2}]⩵0,T[x,0]⩵T0[x],T[0,t]⩵0,T[L,t]⩵0},{T},{x,0,L
Out[133]= T→ InterpolatingFunction Domain: {{0., 10.}, {0., 40.}}
Output: scalar

◼ Získání řešení
In[134]:= fT=T/.dsol[[1]]
Out[134]= InterpolatingFunction Domain: {{0., 10.}, {0., 40.}}
Output: scalar

◼ Zobrazení řešení
symbolicke_vypocty.nb 19
In[135]:= Plot3D[fT[x,t],{x,0,L},{t,0,tmax },AxesLabel→{"x","t","T"}]
Out[135]=
Šíření vln na napnuté membráně
Membrána je tenká deska, která má téměř nulovou tuhost. Nechť membrána s plochou S leží v rovině x-y. Pro příčné kmity
míří výchylka w podél osy z.
Pohybová rovnice vlny šířící se membránou je
∂2
w
∂t2
= c2
∂2
w
∂x2
+
∂2
w
∂y2
. (0.15)
Nechť je membrána trvale v klidu, tj. w=0, tj. funkce w musí splňovat rovnice
w(x, 0, t) = w(a, y, t) = w(x, b, t) = w(0, y, t) = 0. (0.16)
Dále je na počátku pohybu (t=0) známá poloha a rychlost každého bodu membrány, tj.
w(x, y, 0) = f (x, y),
∂w
∂t t=0
= g(x, y), (0.17)
kde f(x,y) a g(x,y) jsou dané funkce místa (souřadnic x,y).
◼ Nechť je rychlost šíření c=1 a rozměry čtvercovémembráby a,b
In[136]:= c=1;
a=10;
b=10;
◼ Iniciální funkce popisující plochu membrány v časet=0 je f(x,y)=Exp(-((x - a)2
+(y - b)2
)). Membrána je uvolněna z klidu, tj.
 dw
dt

0
= g(x, y) = 0.
In[139]:= f[x_,y_]=Exp[-5((x-(a/2)-1)^2+(y-b/2)^2)]+Exp[-1((x-(a/2)+1)^2+(y-b/2)^2)];
20 symbolicke_vypocty.nb
In[140]:= Plot3D[f[x,y],{x,0,a},{y,0,b},PlotRange→Full]
Out[140]=
◼ Numerické řešenívlnové rovnice
In[141]:= tmax =20;
dw[t_,x_,y_]=D[w[t,x,y],t];
dsol=NDSolve[{D[w[t,x,y],{t,2}]-c (D[w[t,x,y],{x,2}]+D[w[t,x,y],{y,2}])⩵0,
w[t,x,0]⩵0,w[t,x,b]⩵0,w[t,0,y]⩵0,w[t,a,y]⩵0,w[0,x,y]⩵f[x,y],
dw[0,x,y]⩵0},{w},{x,0,a},{y,0,b},{t,0,tmax }];
◼ Získané řešení
In[144]:= fw=w/.dsol[[1]]
Out[144]= InterpolatingFunction Domain: {{0., 20.}, {0., 10.}, {0., 10.}}
Output: scalar

◼ Zobrazení tvaru membrány v různých časovýchokamžicích.
In[145]:= t=20;
Plot3D[fw[t,x,y],{x,0,a},{y,0,b},PlotRange→{{0,a},{0,b},{-1,1}}]
Out[146]=
symbolicke_vypocty.nb 21
◼ Animace řešení(Pozor na pomalejších strojích může tato akce trvat pár minut!)
In[147]:=
lst=ListAnimate [Table[Plot3D[fw[t,x,y],{x,0,a},{y,0,b},PlotRange→{{0,a},{0,b},{-1,1}}],{t
In[148]:= lst
Out[148]=
In[224]:= path="/home /kopernik/Lectures/SymbolickeVypocty /Doc/"
Export[path<>"wave.avi",lst,"AVI"]
Out[224]= /home /kopernik/Lectures/SymbolickeVypocty /Doc/
Export::nodir : Directory /home /kopernik /Lectures/SymbolickeVypocty /Doc/ does not exist . 
OpenWrite::noopen : Cannot open /home /kopernik /Lectures/SymbolickeVypocty /Doc/wave.avi . 
Out[225]= $Failed
Řešení Poissonovy rovnice
Poissonova rovnice reprezentuje lokální vyjádření gravitačního zákona a píšeme jej ve tvaru
 ϕ = -4 π ρ. (0.18)
Jestliže gravitační pole neovlivňuje rozložení hmoty, pak nám tato rovnice umožní ze zadaného rozložení hmoty určit výsledné
gravitační pole. Uvažujme následující rozložení hmoty
ρ = 
0.01 pro 0.1 < x < 0.3 & × 0.1 < y < 0.3
0.006 pro - 0.4 < x < -0.2 & - 0.4 < y < -0.2
(0.19)
a určeme výsledný potenciál u(x, y).
22 symbolicke_vypocty.nb
In[149]:= Clear[ρ0];
podmn1 =x>0.1&&x<0.3&&y>0.1&&y<0.3;
podmn2 =x>-0.4&&x<-0.2&&y>-0.4&&y<-0.2;
PoissonEquation=D[u[x,y],{x,2}]+D[u[x,y],{y,2}]⩵-4π ρ[x,y];
ρ[x_,y_]:=Which[podmn1 ,0.01,podmn2 ,0.006,Not[podmn1 ]&&Not[podmn2 ],0];
dsol=NDSolve[{PoissonEquation,u[x,-1]⩵u[x,1]⩵u[-1,y]⩵u[1,y]⩵0},u,{x,-1,1},{y,-1,1}]
Out[154]= u→ InterpolatingFunction Domain: {{-1., 1.}, {-1., 1.}}
Output: scalar

In[155]:= fu=u/.dsol[[1]]
Out[155]= InterpolatingFunction Domain: {{-1., 1.}, {-1., 1.}}
Output: scalar

In[156]:= Plot3D[fu[x,y],{x,-1,1},{y,-1,1},PlotPoints→50]
Out[156]=
Pohyb testovacích částic ve Schwarzchildově metrickém poli
In[234]:=
In[235]:=
symbolicke_vypocty.nb 23