# The hamilton package

[ Tags: bsd3, library, physics, program ] [ Propose Tags ]

## Properties

Versions 0.1.0.0, 0.1.0.1 CHANGELOG.md ad, ansi-wl-pprint, base (>=4.9 && <5), comonad, containers, free, hamilton, hmatrix (>=0.18), hmatrix-gsl (>=0.18), optparse-applicative (>=0.13), typelits-witnesses, vector, vector-sized (>=0.6), vty [details] BSD3 (c) Justin Le 2016 Justin Le justin@jle.im Physics https://github.com/mstksg/hamilton head: git clone https://github.com/mstksg/hamilton Thu Aug 17 08:59:20 UTC 2017 by jle LTSHaskell:0.1.0.0, NixOS:0.1.0.1, Stackage:0.1.0.1, Tumbleweed:0.1.0.0 hamilton-examples 163 total (19 in the last 30 days) (no votes yet) [estimated by rule of succession] λ λ λ Docs available Last success reported on 2017-08-17 Hackage Matrix CI

## Modules

[Index]

#### Maintainer's Corner

For package maintainers and hackage trustees

[back to package description]

# Hamilton

Simulate physics on arbitrary coordinate systems using automatic differentiation and Hamiltonian mechanics. State only an arbitrary parameterization of your system and a potential energy function!

For example, a simulating a double pendulum system by simulating the progression of the angles of each bob:

You only need:

1. 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 half-length
y2 = -cos θ1 - cos θ2 / 2

2. The masses/inertias of each of those cartesian coordinates (m1 for x1 and y1, m2 for x2 and y2)

3. 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 first-order differential equations. This library solves those first-order differential equations for you using automatic differentiation and some matrix manipulation.

See a blog post I wrote on this, and also the hackage documentation and the example runner user guide (and its source).

### 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

-- some helper patterns to pattern match on sized vectors
pattern V2 :: a -> a -> V.Vector 2 a
pattern V2 x y <- (V.toList->[x,y])
where
V2 x y = fromJust (V.fromList [x,y])

pattern V4 :: a -> a -> a -> a -> V.Vector 4 a
pattern V4 x y z a <- (V.toList->[x,y,z,a])
where
V4 x y z a = fromJust (V.fromList [x,y,z,a])


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 configuration-space representation of the state into a phase-space 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 step-by-step:

evolution :: [Phase 2]
evolution = iterate (stepHam 0.1 doublePendulum) phase0


And you can get the position of the coordinates as:

positions :: [R 2]
positions = phsPositions <$> evolution  And the position in the underlying cartesian space as: positions' :: [R 4] positions' = underlyingPos doublePendulum <$> positions


## Example App runner

(Source)

Installation:

$git clone https://github.com/mstksg/hamilton$ cd hamilton
$stack install  Usage: $ hamilton-examples [EXAMPLE] (options)
$hamilton-examples --help$ hamilton-examples [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

### Two-body system under gravity

1. 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 θ

2. The masses/inertias are again m1 for x1 and y1, and m2 for x2 and y2

3. 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 two-body system, assuming m1 is 100 and m2 is 1:

twoBody :: System 4 2
twoBody = mkSystem masses coordinates potential
where
masses :: R 4
masses = vec4 100 100 1 1
coordinates
:: Floating a
=> V.Vector 2 a
-> V.Vector 4 a
coordinates (V2 r θ) = V4 (r1 * cos θ) (r1 * sin θ)
(r2 * cos θ) (r2 * sin θ)
where
r1 =   r *   1 / 101
r2 = - r * 100 / 101
potential
:: Num a
=> V.Vector 4 a
-> a
potential (V2 r _) = - 100 / r


## Potential improvements

• Time-dependent systems: Shouldn't be an problem in theory/math; just add a time parameter before all of the functions. This opens a lot of doors, like deriving inertial forces for free (like the famous Coriolis force and centrifugal force).

The only thing is that it makes the API pretty inconvenient, because it'd require all of the functions to also take a time parameter. Of course, the easy way out/ugly solution would be to just offer two versions of the same function (one for time-independent systems and one for time-dependent systems. But this is un-ideal.

• Velocity-dependent potentials: Would give us the ability to model systems with velocity-dependent Lagrangians like a charged particle in an electromagnetic field, and also dissipative systems, like systems with friction (dependent on signum v) and linear & quadratic wind resistance.

This issue is much harder, theoretically. It involves inverting arbitrary functions forall a. RealFloat a => V.Vector n a -> V.Vector m a. It might be possible with the help of some bidirectionalization techniques, but I can't get the bff package to compile, and I'm not sure how to get bff-mono to work with numeric functions.

If anyone is familiar with bidirectionalization techniques and is willing to help out, please send me a message or open an issue! :)