Automatický výpočet Hamiltonových rovnic z lagranžiánu v pythonu pomocí sympy - vyzkoušené příklady.
+
+
+Tento jednoduchý skript řeší poměrně otravný problém, který
+je však snadno algoritmizovatelný. Jsou to jen opakované parciální derivace,
+jejichž výsledky jsou dosazovány do daných výrazů. Protože se přitom člověk
+snadno splete (hlavně ve znaménku), je lepší to nechat na programu.
+Stačí nadefinovat potřebné symboly pro konstanty a časově závislé proměnné
+(zobecněné souřadnice), pomocí nich pak vyjádřit lagranžián a o zbytek se postará python.
+
+
+
+
\ No newline at end of file
diff --git a/Hamilton/index.html b/Hamilton/index.html
new file mode 100644
index 0000000..3e4bb75
--- /dev/null
+++ b/Hamilton/index.html
@@ -0,0 +1,63 @@
+
+
+
+
+ SOLVER
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Solver pro soustavu Hamiltonových rovnic.
+
Nebudu se tady zabývat tím, co to jsou Hamiltonovy rovnice, na to je teoretická mechanika.
+ Pro tento účel je to soustava diferenciálních rovnic, převedených do tvaru
+ $$ \dot q_1 = f_1 (q_1, p_1, ..., q_i, p_i) $$
+ $$ ... $$
+ $$ \dot q_i = f_i (q_1, p_1, ..., q_i, p_i) $$
+ $$ \dot p_1 = g_1 (q_1, p_1, ..., q_i, p_i) $$
+ $$ ... $$
+ $$ \dot p_i = g_i (q_1, p_1, ..., q_i, p_i) $$
+ kde $ \dot q_i, \dot p_i $ jsou derivace zobecněných souřadnic a hybností podle času. Celkově
+ je to tedy sudý počet rovnic 2*i. Tyto rovnice zapíšeme jednoduše ve správném pořadí a se správným
+ pořadím proměnných ve tvaru, který používá javascript, stanovíme 2*i počátečních podmínek a stanovíme
+ co se má vykreslit. Osy to nemá, každý graf jde od minima do maxima, je to jen pro představu, jak funkce
+ vypadá. Objekt, vstupující do rovnic má tvar {t:t, w:[q1, q2 ... qn, p1, p2, ... pn]}, t značí čas.
+ Default příklad je harmonický oscilátor, lze použít i
+ nelineární oscilátor nebo pro příklad
+ Lorenzův atraktor a nakonec
+ railgun s disipací energie.
+ Stačí kliknout na zvýrazněný text
+ a stisknout tlačítko Submit. Parametry simulace lze jednoduše měnit. Připadá mi to o hodně jednodušší
+ popsat takto problém textem než vymýšlet celé IDE s mnoha parametry. Parametry funkce solve (krok, počet_kroků, počáteční_podmínky, seznam_funkcí),
+ plot (výsledek_ze_solve, popis_os, barva). Objekt popis_os má tvar {x:[list], y:[list]}, přičemž list je buď 't' nebo 'w', číslo_proměnné.
+
+
Vlastní solver používá metodu Runge-Kutta 4.řádu, je poměrně přesná a není to o moc složitější
+ než prostý Euler. Pokud hledáte ty Hamiltonovy rovnice, pomůže tento skript,
+ jeho výstup je pak zhruba tento.
+
#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+description="""
+Tento jednoduchý skript řeší poměrně otravný problém, který
+je však snadno algoritmizovatelný. Jsou to jen opakované parciální derivace,
+jejichž výsledky jsou dosazovány do daných výrazů. Protože se přitom člověk
+snadno splete (hlavně ve znaménku), je lepší to nechat na programu.
+Stačí nadefinovat potřebné symboly pro konstanty a časově závislé proměnné
+(zobecněné souřadnice), pomocí nich pak vyjádřit lagranžián a o zbytek se postará python.
+"""
+fromsympy.coreimportmul,add
+fromsympy.physics.vector.printingimportvlatex
+fromsympy.physics.mechanicsimportdynamicsymbols
+fromsympyimport(Symbol,symbols,Derivative,
+ Eq,pprint,solve,latex,simplify,expand)
+
+################# HTML dekorace ###############################################
+classTag:
+ def__init__(self,n,v):
+ self.name=n
+ self.value=v
+ defto_str(self):
+ return' '+self.name+'=\"'+self.value+'\"'
+
+classElement:
+ def__init__(self,n,v=None):
+ self.name=n
+ self.value=v
+ self.tags=[]
+ self.childs=[]
+ self.id=0
+ defaddE(self,e):
+ self.childs.append(e)
+ defaddT(self,t):
+ self.tags.append(t)
+ defaddF(self,f):
+ t=Tag('class','formulaDsp')
+ forpinf:
+ n=Element('p','\\['+p+'\\]')
+ n.addT(t)
+ self.addE(n)
+
+ defto_str(self):
+ s=self.indent()
+ s+='<'+self.name
+ forninself.tags:s+=n.to_str()
+ s+='>'
+ ifself.value!=None:s+=self.value
+ forninself.childs:s+=n.to_str()
+ s+=self.indent()
+ s+='</'+self.name+'>'
+ returns
+ defcal(self,k=0):
+ self.id=k
+ forninself.childs:n.cal(k+2)
+ defindent(self):
+ s='\n'+self.id*' '
+ returns
+
+defset_root(e,root):
+ styl=Element('style','body {background-color: rgb(192,255,255);} h2 {color: rgb(64,0,192);} h3 {color: rgb(192,0,0);} table {color: rgb(128,0,128);}')
+ scfg='MathJax.Hub.Config({\n extensions: [\"tex2jax.js\","TeX/AMSmath.js"],\n jax: [\"input/TeX\",'
+ scfg+='\"output/HTML-CSS\"],\n tex2jax: {inlineMath: [[\'$\',\'$\']]},\n displayAlign: \"left\"});'
+ head=Element('head')
+ meta=Element('META')
+ titl=Element('title','Lagrange')
+ s1=Element('script',scfg)
+ s2=Element('script')
+ s1.addT(Tag('type','text/x-mathjax-config'))
+ s2.addT(Tag('type','text/javascript'))
+ s2.addT(Tag('src','https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js'))
+ meta.addT(Tag('HTTP-EQUIV','CONTENT-TYPE'))
+ meta.addT(Tag('CONTENT','text/html; charset=utf-8'))
+ head.addE(meta)
+ head.addE(titl)
+ head.addE(styl)
+ head.addE(s1)
+ head.addE(s2)
+ e.addE(head)
+ body=Element('body')
+ body.addE(root)
+ e.addE(body)
+ e.cal()
+defhtml_head():
+ s='<!DOCTYPE html PUBLIC \"-//W3C//DTD XHTML 1.0 Transitional//EN\"\"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd\">'
+ returns
+
+defCreateHTML(filename,html):
+ s=html_head()
+ r=Element('html')
+ set_root(r,html)
+ s+=r.to_str()
+ file=open(filename,'w')
+ file.write(s)
+ file.close()
+
+defoutput(counter,problem,L,H,hce,lce):
+ html=Element('div','<hr>')
+ html.addE(Element('h2','Řešený problém : {0:s}'.format(problem)))
+ html.addE(Element('h3','{0:d}.1. Zadaný lagranžián :'.format(counter)))
+ html.addF(['L = '+vlatex(L)]);
+ html.addE(Element('h3','{0:d}.2. Lagrangeovy rovnice :'.format(counter)))
+ html.addF([vlatex(lce)]);
+ html.addE(Element('h3','{0:d}.3. Hamiltonova funkce :'.format(counter)))
+ html.addF(['H = '+vlatex(H)]);
+ html.addE(Element('h3','{0:d}.4. Hamiltonovy kanonické rovnice :'.format(counter)))
+ table=Element('table')
+ table.addT(Tag('align','center'))
+ forpinhce:
+ td=''
+ forqinp:td+='\n<td class="formulaDsp">\\[{0:s}={1:s}\qquad\\]</td>'.format(vlatex(q[0]),vlatex(expand(q[1])))
+ tr=Element('tr',td)
+ table.addE(tr)
+ html.addE(table)
+ returnhtml
+#################### Vlastní výpočty ##########################################
+t=Symbol('t')# Globální symbol pro čas
+# vykrácení konstantou uděláme ručně, ale nemusí to fungovat, v podstatě jde jen o zbytečné hmotnosti
+defparse(ex):
+ ex=expand(ex)
+ ifex.func!=add.Add:returnex
+ e1=ex.args[0]
+ ife1.func!=mul.Mul:returnex
+ p=e1.args
+ x=[]
+ forbinp:x.append(0)
+ foreinex.args:
+ ife.func!=mul.Mul:break
+ foraine.args:
+ fori,binenumerate(p):
+ ifa==b:x[i]+=1# TODO: or plus power of expression
+ l=len(ex.args)
+ z=[]
+ fori,yinenumerate(x):
+ ify==l:z.append(p[i])
+ foreinz:ex=ex/e
+ ex=simplify(ex)
+ ex=expand(ex)
+ returnex
+
+deflagrange(L,coord):
+ lce=[]
+ forcincoord:
+ ex=L.diff(c[1]).diff(t)-L.diff(c[0])# výraz pro Lagrangeovu rovnici
+ ex=parse(ex)# vykrátit případné konstanty (lze i ručně)
+ le=Eq(ex,0)# a udělat z toho rovnici
+ lce.append(le)
+ pprint(lce)
+ returnlce
+
+defcompute(input,Counter):
+ res=input()
+ coord__x=res[0]
+ coord__p=res[1]
+ L=res[2]
+ pprint(L)# Pro kontrolu
+ coord_dp=[]
+ forpincoord__p:coord_dp.append(p.diff(t))
+ coord_dx=[]
+ forpincoord__x:coord_dx.append(p.diff(t))
+ lce=lagrange(L,zip(coord__x,coord_dx))
+ E=-L
+ forpincoord_dx:E+=L.diff(p)*p
+ E=simplify(E)
+ #pprint (E) # Energie
+ coord_ps=[]
+ forpincoord_dx:coord_ps.append(L.diff(p))
+ coord__h=zip(coord__p,coord_ps)
+ eqs=[]# Legendreova duální transformace
+ forpincoord__h:eqs.append(Eq(p[0],p[1]))
+ sol=[]
+ fori,pinenumerate(eqs):sol.append(solve([p],[coord_dx[i]]))
+ fori,pinenumerate(sol):E=E.subs(coord_dx[i],p[coord_dx[i]])
+ H=simplify(E)
+ H=expand(H)# O něco čitelněji
+ pprint(H)# Hamiltonova funkce
+ hce=[]
+ fori,pinenumerate(coord__x):
+ eqx=Eq(+H.diff(coord__p[i])-coord_dx[i],0)
+ eqx=solve([eqx],[coord_dx[i]])
+ eqx=[coord_dx[i],eqx[coord_dx[i]]]
+ eqp=Eq(-H.diff(coord__x[i])-coord_dp[i],0)
+ eqp=solve([eqp],[coord_dp[i]])
+ eqp=[coord_dp[i],eqp[coord_dp[i]]]
+ hce.append([eqx,eqp])
+ pprint(hce)# Hamiltonovy kanonické rovnice
+ returnoutput(Counter,res[3],L,H,hce,lce)
+
+######################## Uživatelská část #####################################
+defentry1():
+ ##### Zadání #####
+ m,mu,C=symbols('m mu C')# konstanty
+ x,Q=dynamicsymbols('x Q')# proměnné (zobecněné souřadnice)
+ px,pQ=dynamicsymbols('p_x p_Q')# symboly pro příslušné zobecněné hybnosti
+ # Lagranžián
+ L=(m*(x.diff(t))**2)/(2)+(mu*x*(Q.diff(t))**2)/(2)-(Q**2)/(2*C)
+ #############################################################################
+ return[[x,Q],[px,pQ],L,'Railgun']# Zobecníme výstup
+defentry2():# Z učebnice - pohyb planety kolem Slunce, vychází
+ G,m,M=symbols('G m M')# konstanty
+ ##### Zadání #####
+ r,phi=dynamicsymbols('r phi')# proměnné (zobecněné souřadnice)
+ pr,pphi=dynamicsymbols('p_r p_phi')# symboly pro příslušné zobecněné hybnosti
+ # Lagranžián
+ L=(m/2*((r.diff(t))**2+(r**2)*((phi.diff(t))**2)))+((G*m*M)/r)
+ #############################################################################
+ return[[r,phi],[pr,pphi],L,'Pohyb planety']# Zobecníme výstup
+defentry3():# Z učebnice - harmonický oscilátor, vychází
+ ##### Zadání #####
+ m,omega=symbols('m omega')# konstanty
+ x=dynamicsymbols('x')# proměnné (zobecněné souřadnice)
+ px=dynamicsymbols('p_x')# symboly pro příslušné zobecněné hybnosti
+ # Lagranžián
+ L=(m*(x.diff(t))**2-m*omega**2*x**2)/2
+ #############################################################################
+ return[[x],[px],L,'Harmonický oscilátor']# Zobecníme výstup
+###############################################################################
+
+if__name__=="__main__":
+ computed=[entry1,entry2,entry3]# Co vše se bude počítat
+ html=Element('div')
+ html.addE(Element('h1','Automatický výpočet Hamiltonových rovnic z lagranžiánu v pythonu pomocí sympy - vyzkoušené příklady.'))
+ html.addE(Element('p',description))
+ fori,pinenumerate(computed,start=1):html.addE(compute(p,i))
+ CreateHTML('equations.html',html)# převod do lidsky čitelné formy (mathjax latex v html)
+
+
+
+
diff --git a/Hamilton/lorenz.txt b/Hamilton/lorenz.txt
new file mode 100644
index 0000000..396fe46
--- /dev/null
+++ b/Hamilton/lorenz.txt
@@ -0,0 +1,29 @@
+let bc = {t:0.0, w:[0.0, 1.0, 10.0]}
+const omega = 27.0;
+const sigma = 10.0;
+const beta = 8.0 / 3.0;
+let pf = [
+ (q) => {
+ const x = q.w[0];
+ const y = q.w[1];
+ return sigma * (y - x);
+ },
+ (q) => {
+ const x = q.w[0];
+ const y = q.w[1];
+ const z = q.w[2];
+ return x * (omega - z) - y;
+ },
+ (q) => {
+ const x = q.w[0];
+ const y = q.w[1];
+ const z = q.w[2];
+ return x * y - beta * z;
+ },
+];
+let r = solve (0.01, 2000, bc, pf);
+plot (r, {x:['t'], y:['w', 0]}, '#00ff00');
+plot (r, {x:['t'], y:['w', 1]}, '#ff0000');
+plot (r, {x:['t'], y:['w', 2]}, '#0000ff');
+plot (r, {x:['w', 0], y:['w', 2]}, '#00ffff');
+
diff --git a/Hamilton/nvcr.txt b/Hamilton/nvcr.txt
new file mode 100644
index 0000000..24c011e
--- /dev/null
+++ b/Hamilton/nvcr.txt
@@ -0,0 +1,21 @@
+let bc = {t:0.0, w:[0.1, 0.0]}
+const C = 1.0;
+const L = 1.0;
+const G = 0.2;
+const U0= 1.0;
+let pf = [
+ (o) => {
+ const u = o.w[0];
+ const v = o.w[1];
+ return -(v + G*u*(u*u - U0*U0)) / C;
+ },
+ (o) => {
+ const u = o.w[0];
+ return u / L;
+ },
+];
+let r = solve (0.02, 4000, bc, pf);
+plot (r, {x:['t'], y:['w', 0]}, '#00ff00');
+plot (r, {x:['t'], y:['w', 1]}, '#ff0000');
+plot (r, {x:['w', 0], y:['w', 1]}, '#00ffff');
+
diff --git a/Hamilton/osc.txt b/Hamilton/osc.txt
new file mode 100644
index 0000000..909262d
--- /dev/null
+++ b/Hamilton/osc.txt
@@ -0,0 +1,19 @@
+/* Počáteční podmínky */
+let bc = {t:0.0, w:[1.0, 0.0]}
+/* Konstanty */
+const f = 6.2831853
+const fq = f * f;
+/* Hamiltonovy rovnice */
+let pf = [
+ (o) => {
+ return o.w[1];
+ },
+ (o) => {
+ return - fq * o.w[0];
+ },
+];
+let r = solve (0.01, 100, bc, pf); /* řešení */
+/* Grafy */
+plot (r, {x:['t'], y:['w', 0]}, '#00ff00');
+plot (r, {x:['t'], y:['w', 1]}, '#ff0000');
+plot (r, {x:['w', 0], y:['w', 1]}, '#00ffff');
diff --git a/Hamilton/railgun.txt b/Hamilton/railgun.txt
new file mode 100644
index 0000000..2e8b117
--- /dev/null
+++ b/Hamilton/railgun.txt
@@ -0,0 +1,32 @@
+const m = 0.0007, mu = 1.2e-6, C = 0.020, U0 = 250.0, R = 0.01;
+/* {t:t , w:[ x, p_x, Q , p_Q]}*/
+let bc = {t:0.0, w:[0.2, 0.0, C * U0, 0.0]}
+let ls = [
+ (z) => {
+ const p_x = z.w[1];
+ return p_x / m;
+ },
+ (z) => {
+ const x = z.w[0];
+ const p_Q = z.w[3];
+ const d = p_Q / x;
+ return 0.5 * d * d / mu;
+ },
+ (z) => {
+ const x = z.w[0];
+ const p_Q = z.w[3];
+ return p_Q / (mu * x);
+ },
+ (z) => {
+ const x = z.w[0];
+ const Q = z.w[2];
+ const p_Q = z.w[3];
+ const Rp = - (R * p_Q) / (mu * x);
+ return Rp - Q/C;
+ },
+];
+let r = solve (5.0e-7, 2000, bc, ls); /* řešení */
+
+plot (r, {x:['t'], y:['w', 1]}, '#00ff00');
+plot (r, {x:['t'], y:['w', 3]}, '#ff0000');
+
diff --git a/index.html b/index.html
index 0b27c16..48b9986 100644
--- a/index.html
+++ b/index.html
@@ -45,6 +45,10 @@
MIT přikládám. Pro kompilaci je použit jen
clang a jím kompilovaná C-čková knihovna newlib.
+
Že jde něco takového napsat v čistém javascriptu, viz zde.
+ Není to dokonalé, ale funguje to. Python by na to byl asi lepší, ale tady není potřeba nic jiného
+ než prohlížeč.
+