(ns emmy.mechanics.hamilton
(:refer-clojure :exclude [+ - * / partial])
(:require [clojure.core :as core]
[emmy.calculus.derivative :refer [D D-as-matrix partial]]
[emmy.function :as f]
[emmy.generic :as g :refer [sin cos + - * /]]
[emmy.matrix :as matrix]
[emmy.mechanics.lagrange :as l]
[emmy.operator :as o]
[emmy.pattern.rule :as r]
[emmy.structure :as s :refer [up]]
[emmy.value :as v]))

Hamiltonian mechanics requires a phase space QxP, and a function H:RxQxP --> R

A system has a dynamic state, which has the time, the configuration, and the momenta. Hamiltonian mechanics is formulated in terms of the dynamic state.

(defn Hamiltonian
"Returns function signature for a Hamiltonian with n degrees of freedom (or an
unrestricted number if n is not given).
Useful for constructing Hamiltonian literal functions."
[n]
(r/template
(-> (UP Real (UP* Real ~n) (DOWN* Real ~n)) Real)))
#object[emmy.mechanics.hamilton$Hamiltonian 0x1efc23bb "
emmy.mechanics.hamilton$Hamiltonian@1efc23bb"
]
(defn ->H-state
"Given a time `t`, coordinate tuple (or scalar) `q` and momentum tuple (or
scalar) `p`, returns a 'Hamiltonian state tuple', i.e., the state expected by a
Hamiltonian."
[t q p]
(up t q p))
#object[emmy.mechanics.hamilton$__GT_H_state 0x7d7caa3a "
emmy.mechanics.hamilton$__GT_H_state@7d7caa3a"
]
(defn H-state?
"Returns true if the supplied state is
- of type [[emmy.structure/up]]
- contains three elements of `time`, `coordinate` and `momentum` of either of
the following type shapes:
```
(up <number> <number> <number>)
(up <number> (up <number>*) (down <number>*))
```
If structural, the dimension of the coordinate and momentum tuples must match."
[s]
(and (s/up? s)
(= (count s) 3)
(let [[t q v] s]
(and (v/numerical? t)
(or (and (v/numerical? q)
(v/numerical? v))
(and (s/up? q)
(s/down? v)
(= (s/dimension q)
(s/dimension v))))))))
#object[emmy.mechanics.hamilton$H_state_QMARK_ 0x45bfa80c "
emmy.mechanics.hamilton$H_state_QMARK_@45bfa80c"
]
(defn compatible-H-state?
"Returns true if `s` is compatible for contraction with a proper H-state, false
otherwise."
[s]
(H-state?
(s/transpose s)))
#object[emmy.mechanics.hamilton$compatible_H_state_QMARK_ 0x1b3f4af2 "
emmy.mechanics.hamilton$compatible_H_state_QMARK_@1b3f4af2"
]
(defn momentum
"Returns the momentum element of a local Hamiltonian state tuple (by convention,
the third element)."
[H-state]
{:pre [(s/up? H-state) (> (count H-state) 2)]}
(nth H-state 2))
#object[emmy.mechanics.hamilton$momentum 0x642143d5 "
emmy.mechanics.hamilton$momentum@642143d5"
]
(def state->p
"Alias for [[momentum]]."
momentum)
#object[emmy.mechanics.hamilton$momentum 0x642143d5 "
emmy.mechanics.hamilton$momentum@642143d5"
]
(def momenta
"Alias for [[momentum]]."
momentum)
#object[emmy.mechanics.hamilton$momentum 0x642143d5 "
emmy.mechanics.hamilton$momentum@642143d5"
]
(def P
"Alias for [[momentum]]."
momentum)
#object[emmy.mechanics.hamilton$momentum 0x642143d5 "
emmy.mechanics.hamilton$momentum@642143d5"
]
(defn state->qp
"Given a hamiltonian state, returns a [[emmy.structure/up]] containing the
coordinate and momentum components. "
[s]
{:pre [(H-state? s)]}
(up (l/coordinate s)
(momentum s)))
#object[emmy.mechanics.hamilton$state__GT_qp 0x7c141811 "
emmy.mechanics.hamilton$state__GT_qp@7c141811"
]
(defn qp->H-state-path [q p]
(fn [t]
(->H-state t (q t) (p t))))
#object[emmy.mechanics.hamilton$qp__GT_H_state_path 0x44d32057 "
emmy.mechanics.hamilton$qp__GT_H_state_path@44d32057"
]
(defn literal-Hamiltonian-state [n-dof]
(->H-state
(gensym 't)
(s/literal-up (gensym 'x) n-dof)
(s/literal-down (gensym 'p) n-dof)))
#object[emmy.mechanics.hamilton$literal_Hamiltonian_state 0x25da9ae5 "
emmy.mechanics.hamilton$literal_Hamiltonian_state@25da9ae5"
]
(defn L-state->H-state [L]
(fn [Ls]
(->H-state
(l/time Ls)
(l/coordinate Ls)
(((partial 2) L) Ls))))
#object[emmy.mechanics.hamilton$L_state__GT_H_state 0x4b8b4a60 "
emmy.mechanics.hamilton$L_state__GT_H_state@4b8b4a60"
]
(defn H-state->L-state [H]
(fn [Hs]
(l/->L-state
(l/time Hs)
(l/coordinate Hs)
(((partial 2) H) Hs))))
#object[emmy.mechanics.hamilton$H_state__GT_L_state 0x17830790 "
emmy.mechanics.hamilton$H_state__GT_L_state@17830790"
]
(defn H-state->matrix [s]
(matrix/s->m
(s/compatible-shape s) s 1))
#object[emmy.mechanics.hamilton$H_state__GT_matrix 0x4a8f5caf "
emmy.mechanics.hamilton$H_state__GT_matrix@4a8f5caf"
]
(defn matrix->H-state [m s]
{:pre [(= (matrix/num-cols m) 1)
(odd? (matrix/num-rows m))
(> (matrix/num-rows m) 2)]}
(matrix/m->s
(s/compatible-shape s) m 1))
#object[emmy.mechanics.hamilton$matrix__GT_H_state 0x3b73eaa9 "
emmy.mechanics.hamilton$matrix__GT_H_state@3b73eaa9"
]
(defn make-Hamiltonian [kinetic-energy potential-energy]
(+ kinetic-energy potential-energy))
#object[emmy.mechanics.hamilton$make_Hamiltonian 0x2aba6fed "
emmy.mechanics.hamilton$make_Hamiltonian@2aba6fed"
]
(defn Hamiltonian->state-derivative [H]
(fn [H-state]
(->H-state 1
(((partial 2) H) H-state)
(- (((partial 1) H) H-state)))))
#object[emmy.mechanics.hamilton$Hamiltonian__GT_state_derivative 0x258a42ad "
emmy.mechanics.hamilton$Hamiltonian__GT_state_derivative@258a42ad"
]
(def phase-space-derivative
"Alias for [[Hamiltonian->state-derivative]], for compatibility with
1st edition of SICM."
Hamiltonian->state-derivative)
#object[emmy.mechanics.hamilton$Hamiltonian__GT_state_derivative 0x258a42ad "
emmy.mechanics.hamilton$Hamiltonian__GT_state_derivative@258a42ad"
]
(defn Hamilton-equations [Hamiltonian]
(fn [q p]
(let [H-state-path (qp->H-state-path q p)
dH (Hamiltonian->state-derivative Hamiltonian)]
(- (D H-state-path)
(f/compose dH H-state-path)))))
#object[emmy.mechanics.hamilton$Hamilton_equations 0x15d00897 "
emmy.mechanics.hamilton$Hamilton_equations@15d00897"
]
(defn D-phase-space [H]
(fn [s]
(up 0
(((partial 2) H) s)
(- (((partial 1) H) s)))))
#object[emmy.mechanics.hamilton$D_phase_space 0x31d28e26 "
emmy.mechanics.hamilton$D_phase_space@31d28e26"
]

If we express the energy in terms of t,Q,P we have the Hamiltonian. A Hamiltonian is an example of an H-function: an H-function takes 2 vector arguments and a scalar argument (t, Q, P). It produces a scalar result.

(defn H-rectangular [m V]
(fn [[_ q p]]
(make-Hamiltonian
(/ (g/square p) (* 2 m))
(apply V q))))
#object[emmy.mechanics.hamilton$H_rectangular 0x486e4af6 "
emmy.mechanics.hamilton$H_rectangular@486e4af6"
]
(defn H-central [m V]
(fn [[_ q p]]
(make-Hamiltonian
(/ (g/square p) (* 2 m))
(V (g/abs q)))))
#object[emmy.mechanics.hamilton$H_central 0x3f86f67c "
emmy.mechanics.hamilton$H_central@3f86f67c"
]
(defn H-central-polar [m V]
(fn [[_ [r _] [p_r p_phi]]]
(make-Hamiltonian
(/ (+ (g/square p_r)
(g/square (/ p_phi r)))
(* 2 m))
(V r))))
#object[emmy.mechanics.hamilton$H_central_polar 0x2a87d5ab "
emmy.mechanics.hamilton$H_central_polar@2a87d5ab"
]
(defn H-harmonic [m k]
(fn [[_ q p]]
(+ (/ (g/square p) (* 2 m))
(* (/ 1 2) k (g/square q)))))
#object[emmy.mechanics.hamilton$H_harmonic 0xf81dd63 "
emmy.mechanics.hamilton$H_harmonic@f81dd63"
]

If we express the energy in terms of t,Q,P we have the Hamiltonian

H(t,Q,P) = P*Qdot - L(t, Q, Qdot(t, Q, P))

To do this we need to invert P(t, Q, Qdot) to get Qdot(t, Q, P). This is easy when L is a quadratic form in Qdot:

L(t, Q, Qdot) = 1/2*Qdot*M*Qdot + B*Qdot - V

Fortunately this is the case in almost all of Newtonian mechanics, otherwise the P(t,Q,Qdot) function would be much more difficult to invert to obtain Qdot(t,Q,P).

Assume that F is quadratic in its arguments F(u) = 1/2 A u u + b u + c then v = A u + b, so u = A^(-1) (v - b)

(def ^:dynamic *validate-Legendre-transform?*
"If true, the state passed to the fn returned by [[Legendre-transform]] is
checked for correctness. If `false` errors may occur. See the code body for
more detail.
Defaults to `false`."
false)
false
(defn ^:no-doc Legendre-transform-procedure
"Note from GJS: This ugly version tests for correctness of the result."
[F]
(let [w-of-v (D F)
Dw-of-v (D w-of-v)]
(letfn [(putative-G [w]
(let [z (s/compatible-zero w)
M (Dw-of-v z)
b (w-of-v z)]
(if (and *validate-Legendre-transform?*
(g/zero?
(g/simplify
(g/determinant M))))
(do (println "determinant" (g/determinant M))
(throw
(ex-info "Legendre Transform Failure: determinant = 0"
{:F F :w w})))
(let [v (g/solve-linear-left M (- w b))]
(- (* w v) (F v))))))]
(let [Dpg (D putative-G)]
(fn G [w]
(if (and *validate-Legendre-transform?*
(let [thing (s/typical-object w)]
(not (v/= thing
(g/simplify
(w-of-v (Dpg thing)))))))
(throw
(ex-info "Legendre Transform Failure: not quadratic"
{:F F :w w}))
(putative-G w)))))))
#object[emmy.mechanics.hamilton$Legendre_transform_procedure 0x858ea4d "
emmy.mechanics.hamilton$Legendre_transform_procedure@858ea4d"
]
(def Legendre-transform
(o/make-operator Legendre-transform-procedure
'Legendre-transform))
Legendre-transform

Notice that Lagrangians and Hamiltonians are symmetrical with respect to the Legendre transform.

(defn ^:no-doc Lagrangian->Hamiltonian-procedure
[Lagrangian]
(fn [[t q p]] ;; H-state
(let [L #(Lagrangian (up t q %))]
((Legendre-transform L) p))))
#object[emmy.mechanics.hamilton$Lagrangian__GT_Hamiltonian_procedure 0x739c3d75 "
emmy.mechanics.hamilton$Lagrangian__GT_Hamiltonian_procedure@739c3d75"
]
(def Lagrangian->Hamiltonian
(o/make-operator Lagrangian->Hamiltonian-procedure
'Lagrangian->Hamiltonian))
Lagrangian->Hamiltonian
(defn ^:no-doc Hamiltonian->Lagrangian-procedure [Hamiltonian]
(fn [[t q qdot]]
(letfn [(H [p]
(Hamiltonian
(->H-state t q p)))]
((Legendre-transform-procedure H) qdot))))
#object[emmy.mechanics.hamilton$Hamiltonian__GT_Lagrangian_procedure 0x33b25626 "
emmy.mechanics.hamilton$Hamiltonian__GT_Lagrangian_procedure@33b25626"
]
(def Hamiltonian->Lagrangian
(o/make-operator Hamiltonian->Lagrangian-procedure
'Hamiltonian->Lagrangian))
Hamiltonian->Lagrangian
(defn Poisson-bracket [f g]
(fn [x]
(let [fx (f x)
gx (g x)]
(if (or (s/structure? fx) (s/structure? gx))
(s/mapr (fn [af]
(s/mapr (fn [ag]
((Poisson-bracket
(comp (apply s/component af) f)
(comp (apply s/component ag) g))
x))
(s/structure->access-chains gx)))
(s/structure->access-chains fx))
((- (* ((partial 1) f) ((partial 2) g))
(* ((partial 2) f) ((partial 1) g)))
x)))))
#object[emmy.mechanics.hamilton$Poisson_bracket 0x3dbf3a45 "
emmy.mechanics.hamilton$Poisson_bracket@3dbf3a45"
]
(defn ^:no-doc Lie-derivative
"p. 428
We define the Lie derivative of F, as a derivative-like operator, relative to
the given Hamiltonian-like function, H. Generalization and redefinition in
calculus/Lie.scm
"
[H]
(o/make-operator
(fn [F] (Poisson-bracket F H))
(list 'Lie-derivative H)))
#object[emmy.mechanics.hamilton$Lie_derivative 0x29019ad4 "
emmy.mechanics.hamilton$Lie_derivative@29019ad4"
]
(defn flow-derivative
"the flow derivative generalizes the Lie derivative to allow for time dependent
H and F --- computes the 'time' derivative of F along the flow specified by H"
[H]
(-> (fn [F]
(+ ((partial 0) F)
(Poisson-bracket F H)))
(o/make-operator
(list 'flow-derivative H))))
#object[emmy.mechanics.hamilton$flow_derivative 0x286a2052 "
emmy.mechanics.hamilton$flow_derivative@286a2052"
]
(defmethod g/Lie-derivative [::v/function] [f]
(Lie-derivative f))
#object[clojure.lang.MultiFn 0xd31ac27 "
clojure.lang.MultiFn@d31ac27"
]
(defn Lie-transform
"p. 428, the Lie transform is just the time-advance operator using the Lie
derivative (see Hamiltonian.scm)."
[H t]
(-> (g/exp (* t (g/Lie-derivative H)))
(o/make-operator
`(~'Lie-transform ~H ~t))))
#object[emmy.mechanics.hamilton$Lie_transform 0x4998c2d6 "
emmy.mechanics.hamilton$Lie_transform@4998c2d6"
]
(defn flow-transform
"The generalization of Lie-transform to include time dependence."
[H delta-t]
(-> (g/exp (* delta-t (flow-derivative H)))
(o/make-operator
`(~'flow-transform ~H ~delta-t))))
#object[emmy.mechanics.hamilton$flow_transform 0x6cea260b "
emmy.mechanics.hamilton$flow_transform@6cea260b"
]
(defn standard-map [K]
(let [pv (v/principal-value v/twopi)]
(fn [theta I cont _fail]
(let [nI (pv (core/+ I (core/* K (Math/sin theta))))]
(cont
(pv (core/+ theta nI))
nI)))))
#object[emmy.mechanics.hamilton$standard_map 0x157be374 "
emmy.mechanics.hamilton$standard_map@157be374"
]
(defn standard-map-inverse [K]
(let [pv (v/principal-value v/twopi)]
(fn [theta I cont _fail]
(let [ntheta (pv (core/- theta I))]
(cont ntheta
(pv (core/- I (core/*
K (Math/sin ntheta)))))))))
#object[emmy.mechanics.hamilton$standard_map_inverse 0x152b2568 "
emmy.mechanics.hamilton$standard_map_inverse@152b2568"
]
(defn iterated-map
"f is a function of (x y continue fail), which calls continue with the values of
x' y' that follow x y in the mapping.
Returns a map of the same shape that iterates the iterated map n times before
invoking the continuation, or invokes the fail continuation if the inner map
fails."
[f n]
{:pre [(not (neg? n))]}
(let [lulz (constantly nil)]
(fn [x y continue fail]
(loop [x x
y y
i n]
(if (= i 0)
(continue x y)
(if-let [[x' y'] (f x y vector lulz)]
(recur x' y' (dec i))
(fail)))))))
#object[emmy.mechanics.hamilton$iterated_map 0x1a2adb69 "
emmy.mechanics.hamilton$iterated_map@1a2adb69"
]

Point Transformations

Makes a canonical point transformation from a time-invariant coordinate transformation T(q)

(defn F->CH
"A transformation of configuration coordinates F to a procedure implementing a
transformation of phase-space coordinates (p. 320)"
[F]
(fn [[t _ p :as H-state]]
(up t
(F H-state)
(g/solve-linear-right
p
(((partial 1) F) H-state)))))
#object[emmy.mechanics.hamilton$F__GT_CH 0x274d6c59 "
emmy.mechanics.hamilton$F__GT_CH@274d6c59"
]
(def F->CT
"Alias for [[F->CH]]."
F->CH)
#object[emmy.mechanics.hamilton$F__GT_CH 0x274d6c59 "
emmy.mechanics.hamilton$F__GT_CH@274d6c59"
]

This is used in conjunction with a symplectic test for the C to establish that a time-dependent transformation is canonical.

To compute the K (addition to the Hamiltonian) from a time-dependent coordinate transformation F.

(defn F->K [F]
(fn [H-state]
(- (* (g/solve-linear-right
(momentum H-state)
(((partial 1) F) H-state))
(((partial 0) F) H-state)))))
#object[emmy.mechanics.hamilton$F__GT_K 0x6de74b87 "
emmy.mechanics.hamilton$F__GT_K@6de74b87"
]

Canonical Transformations

(defn canonical?
"p.324"
[C H Hprime]
(- (f/compose (Hamiltonian->state-derivative H) C)
(* (D C) (Hamiltonian->state-derivative Hprime))))
#object[emmy.mechanics.hamilton$canonical_QMARK_ 0x21ab3336 "
emmy.mechanics.hamilton$canonical_QMARK_@21ab3336"
]
(defn compositional-canonical?
"p.324"
[C H]
(canonical? C H (f/compose H C)))
#object[emmy.mechanics.hamilton$compositional_canonical_QMARK_ 0x6d9645f9 "
emmy.mechanics.hamilton$compositional_canonical_QMARK_@6d9645f9"
]
(defn J-func [DHs]
(up 0 (nth DHs 2) (- (nth DHs 1))))
#object[emmy.mechanics.hamilton$J_func 0x21f5c3b6 "
emmy.mechanics.hamilton$J_func@21f5c3b6"
]
(defn T-func [s]
(up 1
(g/zero-like (l/coordinates s))
(g/zero-like (momenta s))))
#object[emmy.mechanics.hamilton$T_func 0x76990f7c "
emmy.mechanics.hamilton$T_func@76990f7c"
]
(defn canonical-H? [C H]
(- (f/compose (D-phase-space H) C)
(* (D C)
(D-phase-space (f/compose H C)))))
#object[emmy.mechanics.hamilton$canonical_H_QMARK_ 0xe53a01a "
emmy.mechanics.hamilton$canonical_H_QMARK_@e53a01a"
]
(defn canonical-K? [C K]
(- (f/compose T-func C)
(* (D C)
(+ T-func (D-phase-space K)))))
#object[emmy.mechanics.hamilton$canonical_K_QMARK_ 0x513595a4 "
emmy.mechanics.hamilton$canonical_K_QMARK_@513595a4"
]
(defn linear-function->multiplier [F argument]
((D F) argument))
#object[emmy.mechanics.hamilton$linear_function__GT_multiplier 0x5e521447 "
emmy.mechanics.hamilton$linear_function__GT_multiplier@5e521447"
]
(defn Phi [A]
(fn [v]
(* A v)))
#object[emmy.mechanics.hamilton$Phi 0x608e121e "
emmy.mechanics.hamilton$Phi@608e121e"
]
(defn Phi* [A]
(fn [w]
(* w A)))
#object[emmy.mechanics.hamilton$Phi_STAR_ 0x734e3d65 "
emmy.mechanics.hamilton$Phi_STAR_@734e3d65"
]
(defn time-independent-canonical?
"p.326"
[C]
(fn [s]
((- J-func
(f/compose (Phi ((D C) s))
J-func
(Phi* ((D C) s))))
(s/compatible-shape s))))
#object[emmy.mechanics.hamilton$time_independent_canonical_QMARK_ 0xeab0e1d "
emmy.mechanics.hamilton$time_independent_canonical_QMARK_@eab0e1d"
]

Time-Varying code

Originally from time-varying.scm.

(defn qp-canonical?
"Tests that K yields a canonical transformation if the C is symplectic. (The
qp-canonical? code is really a symplectic test without factoring out the
Hamiltonian.)"
[C H]
(fn [s]
(- (J-func ((D H) (C s)))
(* ((D C) s)
(J-func
((D (f/compose H C)) s))))))
#object[emmy.mechanics.hamilton$qp_canonical_QMARK_ 0x4095c3a3 "
emmy.mechanics.hamilton$qp_canonical_QMARK_@4095c3a3"
]

One particularly useful canonical transform is the Poincare transform, which is good for simplifying oscillators.

(defn polar-canonical
"p.327"
[alpha]
(fn [[t theta I]]
(let [x (* (g/sqrt (/ (* 2 I) alpha)) (sin theta))
p_x (* (g/sqrt (* 2 alpha I)) (cos theta))]
(up t x p_x))))
#object[emmy.mechanics.hamilton$polar_canonical 0x509b374d "
emmy.mechanics.hamilton$polar_canonical@509b374d"
]
(defn polar-canonical-inverse [alpha]
(fn [[t x p]]
(let [I (/ (+ (* alpha (g/square x))
(/ (g/square p) alpha))
2)
theta (g/atan (/ x (g/sqrt (/ (* 2 I) alpha)))
(/ p (g/sqrt (* 2 I alpha))))]
(up t theta I))))
#object[emmy.mechanics.hamilton$polar_canonical_inverse 0x52bc8d2e "
emmy.mechanics.hamilton$polar_canonical_inverse@52bc8d2e"
]
(defn two-particle-center-of-mass [m0 m1]
(fn [[_ [x0 x1]]]
(l/coordinate-tuple
(/ (+ (* m0 x0) (* m1 x1)) (+ m0 m1))
(- x1 x0))))
#object[emmy.mechanics.hamilton$two_particle_center_of_mass 0x4b1aa82 "
emmy.mechanics.hamilton$two_particle_center_of_mass@4b1aa82"
]
(defn two-particle-center-of-mass-canonical [m0 m1]
(fn [[t [x0 x1] [p0 p1]]]
(up t
(l/coordinate-tuple
(/ (+ (* m0 x0) (* m1 x1)) (+ m0 m1))
(- x1 x0))
(l/momentum-tuple
(+ p0 p1)
(/ (- (* m0 p1) (* m1 p0))
(+ m0 m1))))))
#object[emmy.mechanics.hamilton$two_particle_center_of_mass_canonical 0x6ad7f96c "
emmy.mechanics.hamilton$two_particle_center_of_mass_canonical@6ad7f96c"
]
(defn transpose-function [A]
(fn [p]
(* p A)))
#object[emmy.mechanics.hamilton$transpose_function 0x637f9a24 "
emmy.mechanics.hamilton$transpose_function@637f9a24"
]
(defn multiplicative-transpose [s]
(fn [A]
(linear-function->multiplier
(transpose-function A) s)))
#object[emmy.mechanics.hamilton$multiplicative_transpose 0x559cbe98 "
emmy.mechanics.hamilton$multiplicative_transpose@559cbe98"
]

Symplectic

Originally from symplectic.scm.

(defn symplectic-two-form [zeta1 zeta2]
(- (* (momenta zeta2) (l/coordinates zeta1))
(* (momenta zeta1) (l/coordinates zeta2))))
#object[emmy.mechanics.hamilton$symplectic_two_form 0x5cf0eb13 "
emmy.mechanics.hamilton$symplectic_two_form@5cf0eb13"
]

Without matrices

(defn canonical-transform? [C]
(fn [s]
(let [J ((D J-func) (s/compatible-shape s))
DCs ((D C) s)
DCsT (matrix/s:transpose DCs s)]
(- J (* DCs J DCsT)))))
#object[emmy.mechanics.hamilton$canonical_transform_QMARK_ 0x9b97578 "
emmy.mechanics.hamilton$canonical_transform_QMARK_@9b97578"
]
(defn J-matrix
"n == degrees of freedom"
[n]
(let [twon+1 (inc (* 2 n))]
(matrix/generate twon+1 twon+1
(fn [a b]
(cond (zero? a) 0
(zero? b) 0
(= (core/+ a n) b) 1
(= (core/+ b n) a) -1
:else 0)))))
#object[emmy.mechanics.hamilton$J_matrix 0x168870d1 "
emmy.mechanics.hamilton$J_matrix@168870d1"
]
(defn symplectic?
"Symplectic test in terms of matrices"
[C]
(fn [s]
(let [J (J-matrix (l/state->n-dof s))
DC ((D-as-matrix C) s)]
(- J (* DC J (g/transpose DC))))))
#object[emmy.mechanics.hamilton$symplectic_QMARK_ 0x265064e6 "
emmy.mechanics.hamilton$symplectic_QMARK_@265064e6"
]
(defn symplectic-unit
"p. 334 (used, but not defined there)"
[n]
(let [twoN (* 2 n)]
(matrix/generate twoN twoN
(fn [a b]
(cond (= (+ a n) b) 1
(= (+ b n) a) -1
:else 0)))))
#object[emmy.mechanics.hamilton$symplectic_unit 0x7da3ab4c "
emmy.mechanics.hamilton$symplectic_unit@7da3ab4c"
]
(defn symplectic-matrix?
"p. 334"
[M]
(let [two-n (matrix/dimension M)]
(when-not (even? two-n)
(throw
(ex-info "Wrong type -- symplectic-matrix?"
{:M M})))
(let [J (symplectic-unit (quot two-n 2))]
(- J (* M J (g/transpose M))))))
#object[emmy.mechanics.hamilton$symplectic_matrix_QMARK_ 0x3ff31335 "
emmy.mechanics.hamilton$symplectic_matrix_QMARK_@3ff31335"
]
(defn qp-submatrix [m]
(matrix/without m 0 0))
#object[emmy.mechanics.hamilton$qp_submatrix 0x6aa8f2e2 "
emmy.mechanics.hamilton$qp_submatrix@6aa8f2e2"
]
(defn symplectic-transform?
"p. 334"
[C]
(fn [s]
(symplectic-matrix?
(qp-submatrix
((D-as-matrix C) s)))))
#object[emmy.mechanics.hamilton$symplectic_transform_QMARK_ 0x4ffacc05 "
emmy.mechanics.hamilton$symplectic_transform_QMARK_@4ffacc05"
]