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:

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 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 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 = 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:

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

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:

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