hamilton
Physics on generalized coordinate systems using Hamiltonian Mechanics and AD https://github.com/mstksg/hamilton
Version on this page:  [email protected]:2 
LTS Haskell 14.27:  0.1.0.3 
Stackage Nightly 20190921:  0.1.0.3 
Latest on Hackage:  0.1.0.3 
Module documentation for 0.1.0.0
[email protected]:f0a9099cd8474b2fff95c7f027c10b6a3cc43cce1be1a7bebf6914b4a2c2a49d,2104
 Numeric
Hamilton
Simulate physics on arbitrary coordinate systems using automatic differentiation and Hamiltonian mechanics.
For example, a simulating a double pendulum system by simulating the progression of the angles of each bob:
You only need:

Your generalized coordinates (in this case,
θ1
andθ2
), and equations to convert them to cartesian coordinates of your objects:x1 = sin θ1 y1 = cos θ1 x2 = sin θ1 + sin θ2 / 2  second pendulum is halflength y2 = cos θ1  cos θ2 / 2

The masses/inertias of each of those cartesian coordinates (
m1
forx1
andy1
,m2
forx2
andy2
) 
A potential energy function for your objects:
U = (m1 y1 + m2 y2) * g
And that’s it! Hamiltonian mechanics steps your generalized coordinates (θ1
and θ2
) through time, without needing to do any simulation involving
x1
/y1
/x2
/y2
! And you don’t need to worry about tension or any other
stuff like that. All you need is a description of your coordinate system
itself, and the potential energy!
doublePendulum :: System 4 2
doublePendulum =
mkSystem' (vec4 m1 m1 m2 m2)  masses
(\(V2 θ1 θ2) > V4 (sin θ1) (cos θ1)
(sin θ1 + sin θ2/2) (cos θ1  cos θ2/2)
)  coordinates
(\(V4 _ y1 _ y2) > (m1 * y1 + m2 * y2) * g)
 potential
Thanks to ~~Alexander~~ William Rowan Hamilton, we can express our system parameterized by arbitrary coordinates and get back equations of motions as firstorder differential equations. This library solves those firstorder differential equations for you using automatic differentiation and some matrix manipulation.
See documentation and example runner.
Full Exmaple
Let’s turn our double pendulum (with the second pendulum half as long) into an
actual running program. Let’s say that g = 5
, m1 = 1
, and m2 = 2
.
First, the system:
import Numeric.LinearAlgebra.Static
import qualified Data.Vector.Sized as V
doublePendulum :: System 4 2
doublePendulum = mkSystem' masses coordinates potential
where
masses :: R 4
masses = vec4 1 1 2 2
coordinates
:: Floating a
=> V.Vector 2 a
> V.Vector 4 a
coordinates (V2 θ1 θ2) = V4 (sin θ1) (cos θ1)
(sin θ1 + sin θ2/2) (cos θ1  cos θ2/2)
potential
:: Num a
=> V.Vector 4 a
> a
potential (V4 _ y1 _ y2) = (y1 + 2 * y2) * 5
Neat! Easy, right?
Okay, now let’s run it. Let’s pick a starting configuration (state of the
system) of θ1
and θ2
:
config0 :: Config 2
config0 = Cfg (vec2 1 0 )  initial positions
(vec2 0 0.5)  initial velocities
Configurations are nice, but Hamiltonian dynamics is all about motion through phase space, so let’s convert this configurationspace representation of the state into a phasespace representation of the state:
phase0 :: Phase 2
phase0 = toPhase doublePendulum config0
And now we can ask for the state of our system at any amount of points in time!
ghci> evolveHam doublePendulum phase0 [0,0.1 .. 1]
 result: state of the system at times 0, 0.1, 0.2, 0.3 ... etc.
Or, if you want to run the system stepbystep:
evolution :: [Phase 2]
evolution = iterate (stepHam 0.1 doublePendulum) phase0
And you can get the position of the coordinates as:
positions :: [R 2]
positions = phsPos <$> evolution
And the position in the underlying cartesian space as:
positions' :: [R 4]
positions' = underlyingPos doublePendulum <$> positions
Example App runner
Installation:
$ git clone https://github.com/mstksg/hamilton
$ cd hamilton
$ stack install
Usage:
$ hamiltonexamples [EXAMPLE] (options)
$ hamiltonexamples help
$ hamiltonexamples [EXAMPLE] help
The example runner is a command line application that plots the progression of several example system through time.
Example  Description  Coordinates  Options 

doublepend 
Double pendulum, described above  θ1 , θ2 (angles of bobs) 
Masses of each bob 
pend 
Single pendulum  θ (angle of bob) 
Initial angle and velocity of bob 
room 
Object bounding around walled room  x , y 
Initial launch angle of object 
twobody 
Two gravitationally attracted bodies, described below  r , θ (distance between bodies, angle of rotation) 
Masses of bodies and initial angular veocity 
spring 
Spring hanging from a block on a rail, holding up a weight  r , x , θ (position of block, spring compression, spring angle) 
Masses of block, weight, spring constant, initial compression 
bezier 
Bead sliding at constant velocity along bezier curve  t (Bezier time parameter) 
Control points for arbitrary bezier curve 
Call with help
(or [EXAMPLE] help
) for more information.
More examples
Twobody system under gravity

The generalized coordinates are just:
r
, the distance between the two bodiesθ
, the current angle of rotation
x1 = m2/(m1+m2) * r * sin θ  assuming (0,0) is the center of mass y1 = m2/(m1+m2) * r * cos θ x2 = m1/(m1+m2) * r * sin θ y2 = m1/(m1+m2) * r * cos θ

The masses/inertias are again
m1
forx1
andy1
, andm2
forx2
andy2

The potential energy function is the classic gravitational potential:
U =  m1 * m2 / r
And…that’s all you need!
Here is the actual code for the twobody system:
twoBody :: System 4 2
twoBody =
mkSystem (vec4 m1 m1 m2 m2)  masses
(\(V2 r θ) > let r1 = r * m2 / (m1 + m2)
r2 =  r * m1 / (m1 + m2)
in V4 (r1 * cos θ) (r1 * sin θ)
(r2 * cos θ) (r2 * sin θ)
)  coordinates
(\(V2 r _) >  m1 * m2 / r)  potential