diff --git a/Hamilton/equations.html b/Hamilton/equations.html new file mode 100644 index 0000000..6a29198 --- /dev/null +++ b/Hamilton/equations.html @@ -0,0 +1,114 @@ + + + + + + Lagrange + + + + + + +
+

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. + +

+
+

Řešený problém : Railgun +

+

1.1. Zadaný lagranžián : +

+

\[L = \frac{m \left(\dot{x}\right)^{2}}{2} + \frac{\mu \left(\dot{Q}\right)^{2}}{2} x - \frac{Q^{2}}{2 C}\] +

+

1.2. Lagrangeovy rovnice : +

+

\[\left [ m \ddot{x} - \frac{\mu \left(\dot{Q}\right)^{2}}{2} = 0, \quad \mu x \ddot{Q} + \mu \dot{Q} \dot{x} + \frac{Q}{C} = 0\right ]\] +

+

1.3. Hamiltonova funkce : +

+

\[H = \frac{p^{2}_{Q}}{2 \mu x} + \frac{p^{2}_{x}}{2 m} + \frac{Q^{2}}{2 C}\] +

+

1.4. Hamiltonovy kanonické rovnice : +

+ + + + + + + + + +
\[\dot{x}=\frac{p_{x}}{m}\qquad\]\[\dot{p}_{x}=\frac{p^{2}_{Q}}{2 \mu x^{2}}\qquad\]
\[\dot{Q}=\frac{p_{Q}}{\mu x}\qquad\]\[\dot{p}_{Q}=- \frac{Q}{C}\qquad\]
+
+
+

Řešený problém : Pohyb planety +

+

2.1. Zadaný lagranžián : +

+

\[L = \frac{G M}{r} m + \frac{m}{2} \left(r^{2} \left(\dot{\phi}\right)^{2} + \left(\dot{r}\right)^{2}\right)\] +

+

2.2. Lagrangeovy rovnice : +

+

\[\left [ \frac{G M}{r^{2}} - r \left(\dot{\phi}\right)^{2} + \ddot{r} = 0, \quad r^{2} \ddot{\phi} + 2 r \dot{\phi} \dot{r} = 0\right ]\] +

+

2.3. Hamiltonova funkce : +

+

\[H = - \frac{G M}{r} m + \frac{p^{2}_{\phi}}{2 m r^{2}} + \frac{p^{2}_{r}}{2 m}\] +

+

2.4. Hamiltonovy kanonické rovnice : +

+ + + + + + + + + +
\[\dot{r}=\frac{p_{r}}{m}\qquad\]\[\dot{p}_{r}=- \frac{G M}{r^{2}} m + \frac{p^{2}_{\phi}}{m r^{3}}\qquad\]
\[\dot{\phi}=\frac{p_{\phi}}{m r^{2}}\qquad\]\[\dot{p}_{\phi}=0\qquad\]
+
+
+

Řešený problém : Harmonický oscilátor +

+

3.1. Zadaný lagranžián : +

+

\[L = - \frac{m \omega^{2}}{2} x^{2} + \frac{m \left(\dot{x}\right)^{2}}{2}\] +

+

3.2. Lagrangeovy rovnice : +

+

\[\left [ \omega^{2} x + \ddot{x} = 0\right ]\] +

+

3.3. Hamiltonova funkce : +

+

\[H = \frac{m \omega^{2}}{2} x^{2} + \frac{p^{2}_{x}}{2 m}\] +

+

3.4. Hamiltonovy kanonické rovnice : +

+ + + + + +
\[\dot{x}=\frac{p_{x}}{m}\qquad\]\[\dot{p}_{x}=- m \omega^{2} x\qquad\]
+
+
+ + \ 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. +

+
+ + diff --git a/Hamilton/index.js b/Hamilton/index.js new file mode 100644 index 0000000..c67e50d --- /dev/null +++ b/Hamilton/index.js @@ -0,0 +1,126 @@ +var Outs, Inps, Btnt, Canvas; +// async/await z příkladu na webu +window.onload = async function() { + Outs = document.getElementById('stdout'); + Inps = document.getElementById('input'); + Btnt = document.getElementById('buttonTest'); + Canvas = document.getElementById('canvas'); + Outs.value = ''; + await getFile ('osc.txt'); + Btnt.onclick = () => { + Outs.value = ''; + run (); + }; +}; +function print (str) { + Outs.value += str; +} +async function getFile (name) { + console.log(name); + const response = await fetch (name); + if (!response.ok) return; + const bytes = await response.arrayBuffer(); + const array = new Uint8Array(bytes); + const decoder = new TextDecoder(); + Inps.value = decoder.decode(bytes); + return false; +} +/////////////////////////////////////////////// +function run () { + + let code = Inps.value; + try { + eval (code); + } catch (err) { + print('Error: ' + err.message + '\n'); + } +} +function solve (hh, steps, data, pf) { + let p = data; + + const N = pf.length; + const hal = 0.5 * hh; + const hhs = hh / 6.0; + let min = {t:+1.0e100, w:[]}; + for (let k=0; k max.t) max.t = p.t; + for (let k=0; k max.w[k]) max.w[k] = p.w[k]; } + result [0] = new Object (structuredClone(p)); + for (let index=1; index max.t) max.t = np.t; + for (let k=0; k max.w[k]) max.w[k] = np.w[k]; } + result [index] = new Object (np); + } + initCanvas () + return {data:result, m:min, p:max}; +} +function getValue (pt, si) { + if (si.length > 1) { // w + return pt[si[0]][si[1]]; + } else { // t + return pt[si[0]]; + } +} +function initCanvas () { + const width = canvas.clientWidth; + const height = canvas.clientHeight; + canvas.width = width; + canvas.height = height; + const ctx = canvas.getContext("2d"); + ctx.clearRect(0, 0, canvas.width, canvas.height); +} +function plot (r, po, color) { + print('plot:' + JSON.stringify(po) + '\n'); + const p0 = r.data[0]; + print('min: ' + JSON.stringify(r.m) + '\n'); + print('max: ' + JSON.stringify(r.p) + '\n'); + const xmin = getValue (r.m, po.x); + const xmax = getValue (r.p, po.x); + const ymin = getValue (r.m, po.y); + const ymax = getValue (r.p, po.y); + const xzom = canvas.width / (xmax - xmin); + const yzom = canvas.height / (ymin - ymax); + const xofs = xmin * canvas.width / (xmin - xmax); + const yofs = ymax * canvas.height / (ymax - ymin); + //console.log (xmin, xmax, ymin, ymax); + const ctx = canvas.getContext("2d"); + ctx.beginPath(); + const x0 = xzom * getValue(p0, po.x) + xofs; + const y0 = yzom * getValue(p0, po.y) + yofs; + ctx.moveTo(x0, y0); + for (let n=1; n + + + + + + + + +

+ +
  1
+  2
+  3
+  4
+  5
+  6
+  7
+  8
+  9
+ 10
+ 11
+ 12
+ 13
+ 14
+ 15
+ 16
+ 17
+ 18
+ 19
+ 20
+ 21
+ 22
+ 23
+ 24
+ 25
+ 26
+ 27
+ 28
+ 29
+ 30
+ 31
+ 32
+ 33
+ 34
+ 35
+ 36
+ 37
+ 38
+ 39
+ 40
+ 41
+ 42
+ 43
+ 44
+ 45
+ 46
+ 47
+ 48
+ 49
+ 50
+ 51
+ 52
+ 53
+ 54
+ 55
+ 56
+ 57
+ 58
+ 59
+ 60
+ 61
+ 62
+ 63
+ 64
+ 65
+ 66
+ 67
+ 68
+ 69
+ 70
+ 71
+ 72
+ 73
+ 74
+ 75
+ 76
+ 77
+ 78
+ 79
+ 80
+ 81
+ 82
+ 83
+ 84
+ 85
+ 86
+ 87
+ 88
+ 89
+ 90
+ 91
+ 92
+ 93
+ 94
+ 95
+ 96
+ 97
+ 98
+ 99
+100
+101
+102
+103
+104
+105
+106
+107
+108
+109
+110
+111
+112
+113
+114
+115
+116
+117
+118
+119
+120
+121
+122
+123
+124
+125
+126
+127
+128
+129
+130
+131
+132
+133
+134
+135
+136
+137
+138
+139
+140
+141
+142
+143
+144
+145
+146
+147
+148
+149
+150
+151
+152
+153
+154
+155
+156
+157
+158
+159
+160
+161
+162
+163
+164
+165
+166
+167
+168
+169
+170
+171
+172
+173
+174
+175
+176
+177
+178
+179
+180
+181
+182
+183
+184
+185
+186
+187
+188
+189
+190
+191
+192
+193
+194
+195
+196
+197
+198
+199
+200
+201
+202
+203
+204
+205
+206
+207
+208
+209
+210
+211
+212
+213
+214
+215
+216
+217
+218
+219
+220
+221
+222
+223
+224
+225
+226
+227
#!/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.
+"""
+from sympy.core import mul, add
+from sympy.physics.vector.printing import vlatex
+from sympy.physics.mechanics import dynamicsymbols
+from sympy import (Symbol, symbols, Derivative,
+                  Eq, pprint, solve, latex, simplify, expand)
+
+################# HTML dekorace ###############################################
+class Tag:
+  def __init__(self, n, v):
+      self.name  = n
+      self.value = v
+  def to_str (self):
+    return ' ' +  self.name + '=\"' + self.value + '\"'
+
+class Element:
+  def __init__(self, n, v=None):
+      self.name   = n
+      self.value  = v
+      self.tags   = []
+      self.childs = []
+      self.id = 0
+  def addE (self, e):
+    self.childs.append(e)
+  def addT (self, t):
+    self.tags.append(t)
+  def addF (self, f):
+    t = Tag ('class','formulaDsp')
+    for p in f:
+      n = Element('p', '\\[' + p + '\\]')
+      n.addT(t)
+      self.addE(n)
+    
+  def to_str (self):
+    s = self.indent()
+    s+= '<' + self.name
+    for n in self.tags: s += n.to_str()
+    s+= '>'
+    if self.value != None: s += self.value
+    for n in self.childs:  s += n.to_str()
+    s+= self.indent()
+    s+= '</' + self.name + '>'
+    return s
+  def cal (self, k=0):
+    self.id = k
+    for n in self.childs: n.cal (k+2)
+  def indent (self):
+    s = '\n' + self.id * ' '
+    return s
+
+def set_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()
+def html_head ():
+  s = '<!DOCTYPE html PUBLIC \"-//W3C//DTD XHTML 1.0 Transitional//EN\" \"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd\">'
+  return s
+
+def CreateHTML (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()
+
+def output (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'))
+  for p in hce:
+    td = ''
+    for q in p: 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)
+  return html
+#################### 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
+def parse (ex):
+  ex = expand   (ex)
+  if ex.func != add.Add: return ex
+  e1 = ex.args[0]
+  if e1.func != mul.Mul: return ex
+  p = e1.args
+  x = []
+  for b in p: x.append(0)
+  for e in ex.args:
+    if e.func != mul.Mul: break
+    for a in e.args:
+      for i,b in enumerate(p):
+        if a == b: x[i] += 1  # TODO: or plus power of expression
+  l = len (ex.args)
+  z = []
+  for i,y in enumerate(x):
+    if y == l: z.append (p[i])
+  for e in z: ex = ex / e
+  ex = simplify (ex)
+  ex = expand   (ex)
+  return ex
+
+def lagrange (L, coord):
+  lce = []
+  for c in coord:
+    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)
+  return lce
+
+def compute (input, Counter):
+  res = input()
+  coord__x = res[0]
+  coord__p = res[1]
+  L        = res[2]
+  pprint (L)          # Pro kontrolu
+  coord_dp = []
+  for p in coord__p: coord_dp.append (p.diff(t))
+  coord_dx = []
+  for p in coord__x: coord_dx.append (p.diff(t))
+  lce = lagrange (L, zip(coord__x, coord_dx))
+  E = -L
+  for p in coord_dx: E += L.diff (p) * p
+  E  = simplify (E)
+  #pprint (E)         # Energie
+  coord_ps = []
+  for p in coord_dx: coord_ps.append (L.diff(p))
+  coord__h = zip (coord__p, coord_ps)
+  eqs = []            # Legendreova duální transformace
+  for p in coord__h: eqs.append(Eq(p[0], p[1]))
+  sol = []
+  for i,p in enumerate(eqs): sol.append (solve ([p],   [coord_dx[i]]))
+  for i,p in enumerate(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 = []
+  for i,p in enumerate(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
+  return output (Counter, res[3], L, H, hce, lce)
+  
+######################## Uživatelská část #####################################
+def entry1 ():
+  ##### 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
+def entry2 ():   # 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
+def entry3 ():   # 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))
+  for i,p in enumerate(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č. +