\documentclass[]{article}
\usepackage{lmodern}
\usepackage{amssymb,amsmath}
\usepackage{ifxetex,ifluatex}
\usepackage{fixltx2e} % provides \textsubscript
\ifnum 0\ifxetex 1\fi\ifluatex 1\fi=0 % if pdftex
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\else % if luatex or xelatex
\ifxetex
\usepackage{mathspec}
\usepackage{xltxtra,xunicode}
\else
\usepackage{fontspec}
\fi
\defaultfontfeatures{Mapping=tex-text,Scale=MatchLowercase}
\newcommand{\euro}{€}
\fi
% use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
% use microtype if available
\IfFileExists{microtype.sty}{\usepackage{microtype}}{}
\usepackage[margin=1in]{geometry}
\usepackage{color}
\usepackage{fancyvrb}
\newcommand{\VerbBar}{|}
\newcommand{\VERB}{\Verb[commandchars=\\\{\}]}
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}}
% Add ',fontsize=\small' for more characters per line
\newenvironment{Shaded}{}{}
\newcommand{\AlertTok}[1]{\textcolor[rgb]{1.00,0.00,0.00}{\textbf{#1}}}
\newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{#1}}}}
\newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.49,0.56,0.16}{#1}}
\newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{#1}}
\newcommand{\BuiltInTok}[1]{#1}
\newcommand{\CharTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{#1}}
\newcommand{\CommentTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textit{#1}}}
\newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{#1}}}}
\newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.53,0.00,0.00}{#1}}
\newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{\textbf{#1}}}
\newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.56,0.13,0.00}{#1}}
\newcommand{\DecValTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{#1}}
\newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.73,0.13,0.13}{\textit{#1}}}
\newcommand{\ErrorTok}[1]{\textcolor[rgb]{1.00,0.00,0.00}{\textbf{#1}}}
\newcommand{\ExtensionTok}[1]{#1}
\newcommand{\FloatTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{#1}}
\newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.02,0.16,0.49}{#1}}
\newcommand{\ImportTok}[1]{#1}
\newcommand{\InformationTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{#1}}}}
\newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{\textbf{#1}}}
\newcommand{\NormalTok}[1]{#1}
\newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.40,0.40,0.40}{#1}}
\newcommand{\OtherTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{#1}}
\newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.74,0.48,0.00}{#1}}
\newcommand{\RegionMarkerTok}[1]{#1}
\newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{#1}}
\newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.73,0.40,0.53}{#1}}
\newcommand{\StringTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{#1}}
\newcommand{\VariableTok}[1]{\textcolor[rgb]{0.10,0.09,0.49}{#1}}
\newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{#1}}
\newcommand{\WarningTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{#1}}}}
\usepackage{graphicx}
\makeatletter
\def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth\else\Gin@nat@width\fi}
\def\maxheight{\ifdim\Gin@nat@height>\textheight\textheight\else\Gin@nat@height\fi}
\makeatother
% Scale images if necessary, so that they will not overflow the page
% margins by default, and it is still possible to overwrite the defaults
% using explicit options in \includegraphics[width, height, ...]{}
\setkeys{Gin}{width=\maxwidth,height=\maxheight,keepaspectratio}
\ifxetex
\usepackage[setpagesize=false, % page size defined by xetex
unicode=false, % unicode breaks when used with xetex
xetex]{hyperref}
\else
\usepackage[unicode=true]{hyperref}
\fi
\hypersetup{breaklinks=true,
bookmarks=true,
pdfauthor={Justin Le},
pdftitle={Introducing the Hamilton library},
colorlinks=true,
citecolor=blue,
urlcolor=blue,
linkcolor=magenta,
pdfborder={0 0 0}}
\urlstyle{same} % don't use monospace font for urls
% Make links footnotes instead of hotlinks:
\renewcommand{\href}[2]{#2\footnote{\url{#1}}}
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}
\setlength{\emergencystretch}{3em} % prevent overfull lines
\setcounter{secnumdepth}{0}
\title{Introducing the Hamilton library}
\author{Justin Le}
\date{November 28, 2016}
\begin{document}
\maketitle
\emph{Originally posted on
\textbf{\href{https://blog.jle.im/entry/introducing-the-hamilton-library.html}{in
Code}}.}
\href{http://i.imgur.com/Vaaa2EC.gifv}{\includegraphics{/img/entries/hamilton/double-pendulum.gif}}
\textbf{Hamilton}: \href{https://github.com/mstksg/hamilton\#readme}{README} /
\href{http://hackage.haskell.org/package/hamilton}{hackage} /
\href{https://github.com/mstksg/hamilton}{github}
The \emph{\href{http://hackage.haskell.org/package/hamilton}{hamilton}} library
is on hackage! It was mostly a proof-of-concept toy experiment to simulate
motion on bezier curves, but it became usable enough and accurate enough (to my
surprise, admittedly) that I finished up some final touches to make it complete
and put it on hackage as a general-purpose physics simulator.
The library is, in short, a way to simulate a physical system by stating nothing
more than an arbitrary parameterization of a system (a ``generalized
coordinate'') and a potential energy function.
I was going to write a Haskell post on the implementation, which was what
interested me at first. I wanted to go over --
\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
Using \href{https://hackage.haskell.org/package/ad}{automatic differentiation}
to automatically compute momentum and the hamilton equations, which are
solutions of differential equations.
\item
Using type-indexed vectors and dependent types in a seamless way to encode the
dimensionality of the generalized coordinate systems and to encode invariants
the types of functions.
\end{enumerate}
And fun stuff like that. But that post might be a bit of a while away, so I'm
just going to write a post about the usage of the library. (Fair warning, most
of this information is also found in the
\href{https://github.com/mstksg/hamilton\#readme}{readme}.)
\hypertarget{hamiltonian-mechanics}{%
\subsection{Hamiltonian Mechanics}\label{hamiltonian-mechanics}}
\emph{(This section goes briefly over some relevant part of the physics behind
Hamiltonian dynamics, but feel free to skip it if you want to go straight to the
Haskell)}
\href{https://en.wikipedia.org/wiki/Hamiltonian_mechanics}{Hamiltonian
mechanics} is a brilliant, radical, and beautiful re-imagination of the physics
of mechanics and dynamics by
\href{https://www.youtube.com/watch?v=SZXHoWwBcDc}{William Rowan Hamilton}. It
was adapted for statistical mechanics and thermodynamics, and it was through the
lens of Hamiltonian mechanics that Schroedinger and Heisenberg independently
found insight that unlocked the secrets of quantum mechanics. While Newton's
interpretation of mechanics (in terms of forces and accelerations) was cute, it
simply didn't generalize to quantum mechanics. Hamiltonian's interpretation of
mechanics \emph{did}, and we have a century of physics revolutions to thank for
it. Hamiltonian mechanics also generalize without any extra work to relativity
-- another case where newtonian mechanics tends to fall apart.
Hamiltonian mechanics, in a classical sense, imagines that the state of the
system exists as a point in
\emph{\href{https://en.wikipedia.org/wiki/Phase_space}{phase space}}, and that
the system evolves based on geometric properties of the system's
\emph{Hamiltonian} over that phase space.
\begin{figure}
\centering
\includegraphics{/img/entries/hamilton/phase-space.gif}
\caption{Animation of particles traveling in phase space (top) over time, from
Wikipedia}
\end{figure}
In other words, define the Hamiltonian of the system, and you see the
step-by-step evolution and dynamics of the system. You can imagine mechanics as
a series of streams of flow over phase space\ldots and the state of the system
just goes along for the ride.
One nice thing about phase space is that it can be stated in terms of any
arbitrary parameterization/coordinate system of your system. For example, for a
\href{https://en.wikipedia.org/wiki/Double_pendulum}{double pendulum} system,
you can imagine the system as traveling about in the phase space of the angles
of the bobs (instead of their actual positions in cartesian space). If you can
find \emph{any} way to parameterize your system, in any sort of type of
coordinates, then Hamiltonian mechanics will describe how it evolves in those
coordinates.
State some fundamental geometric properties about your coordinate system, and
the Hamiltonian figures out the rest. It's the key to unlocking the dynamical
properties of the system.
I could go into more details, but this isn't a post about Hamiltonian mechanics!
Armed with this, let's look into modeling an actual double pendulum system in
terms of the angles of the bobs.
\hypertarget{examples}{%
\subsection{Examples}\label{examples}}
\hypertarget{the-double-pendulum}{%
\subsubsection{The Double Pendulum}\label{the-double-pendulum}}
So, if we're going to be simulating a double pendulum system using
\emph{hamilton}, we need three things:
\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
A statement of our parameterized coordinates and how they relate to the
underlying cartesian coordinates of our system
\item
The inertias (in our case, masses) of each of those underlying coordinates.
\item
A potential energy function (in our case, just the potential energy induced by
gravity)
\end{enumerate}
We have two coordinates here (\(\theta_1\) and \(\theta_2\)), which will be
encoding the positions of the two pendulums:
\[
\langle x_1, y_1 \rangle =
\left\langle \sin (\theta_1), - \cos (\theta_1) \right\rangle
\]
\[
\langle x_2, y_2 \rangle =
\left\langle \sin (\theta_1) + \frac{1}{2} \sin (\theta_2),
- \cos (\theta_1) - \frac{1}{2} \cos (\theta_2) \right\rangle
\]
(Assuming that the first pendulum has length 1 and the second pendulum has
length \(\frac{1}{2}\))
The inertias of \(x_1\), \(y_1\), \(x_2\), and \(y_2\) are the ``masses''
attached to them. Let's pick that the first bob has mass \(1\) and the second
bob has mass \(2\), so then our masses are \(\langle 1, 1, 2, 2 \rangle\).
Finally, the potential energy of our system is just the potential energy of
gravity, \(m \times g \times y\) for each of our points:
\[
U(x_1, y_1, x_2, y_2) = ( y_1 + 2 y_2 ) g
\]
Turns out that this is a complete enough description of our system to let
\emph{hamilton} do the rest!
\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{{-}{-} source: https://github.com/mstksg/inCode/tree/master/code{-}samples/hamilton/DoublePendulum.hs\#L10{-}L25}
\OtherTok{doublePendulum ::} \DataTypeTok{System} \DecValTok{4} \DecValTok{2}
\NormalTok{doublePendulum }\OtherTok{=}\NormalTok{ mkSystem\textquotesingle{} masses coordinates potential}
\KeywordTok{where}
\OtherTok{ masses ::} \DataTypeTok{R} \DecValTok{4}
\NormalTok{ masses }\OtherTok{=}\NormalTok{ vec4 }\DecValTok{1} \DecValTok{1} \DecValTok{2} \DecValTok{2}
\NormalTok{ coordinates}
\OtherTok{ ::} \DataTypeTok{Floating}\NormalTok{ a}
\OtherTok{=>} \DataTypeTok{V.Vector} \DecValTok{2}\NormalTok{ a}
\OtherTok{{-}>} \DataTypeTok{V.Vector} \DecValTok{4}\NormalTok{ a}
\NormalTok{ coordinates (}\DataTypeTok{V2}\NormalTok{ θ1 θ2) }\OtherTok{=} \DataTypeTok{V4}\NormalTok{ (}\FunctionTok{sin}\NormalTok{ θ1) (}\OperatorTok{{-}}\FunctionTok{cos}\NormalTok{ θ1)}
\NormalTok{ (}\FunctionTok{sin}\NormalTok{ θ1 }\OperatorTok{+} \FunctionTok{sin}\NormalTok{ θ2}\OperatorTok{/}\DecValTok{2}\NormalTok{) (}\OperatorTok{{-}}\FunctionTok{cos}\NormalTok{ θ1 }\OperatorTok{{-}} \FunctionTok{cos}\NormalTok{ θ2}\OperatorTok{/}\DecValTok{2}\NormalTok{)}
\NormalTok{ potential}
\OtherTok{ ::} \DataTypeTok{Num}\NormalTok{ a}
\OtherTok{=>} \DataTypeTok{V.Vector} \DecValTok{4}\NormalTok{ a}
\OtherTok{{-}>}\NormalTok{ a}
\NormalTok{ potential (}\DataTypeTok{V4}\NormalTok{ \_ y1 \_ y2) }\OtherTok{=}\NormalTok{ (y1 }\OperatorTok{+} \DecValTok{2} \OperatorTok{*}\NormalTok{ y2) }\OperatorTok{*} \DecValTok{5} \CommentTok{{-}{-} assuming g = 5}
\end{Highlighting}
\end{Shaded}
(with some
\href{https://github.com/mstksg/inCode/tree/master/code-samples/hamilton/DoublePendulum.hs\#L27-L35}{helper
patterns} defined here -- \texttt{V2} and \texttt{V4} -- that lets us pattern
match on and construct sized \texttt{Vector}s and their 2 (or 4) elements)
Ta dah. That's literally all we need.
A \texttt{System\ m\ n} represents a description of a physical system (without
its state) described with \texttt{n} parameters/generalized coordinates. The
\texttt{m} represents the dimension of its underlying cartesian coordinate
system (\texttt{4} for us, with \(\langle x_1, y_1, x_2, y_2 \rangle\)). The
\texttt{m} should be more or less irrelevant to the actual \emph{usage} of
\texttt{System\ m\ n} and the \emph{hamilton} api\ldots but it's mostly useful
only if we eventually want to plot the system in normal cartesian space.
Now, let's run the simulation. First we have to pick a starting configuration:
\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{{-}{-} source: https://github.com/mstksg/inCode/tree/master/code{-}samples/hamilton/DoublePendulum.hs\#L37{-}L39}
\OtherTok{config0 ::} \DataTypeTok{Config} \DecValTok{2}
\NormalTok{config0 }\OtherTok{=} \DataTypeTok{Cfg}\NormalTok{ (vec2 }\DecValTok{1} \DecValTok{0}\NormalTok{ ) }\CommentTok{{-}{-} initial positions}
\NormalTok{ (vec2 }\DecValTok{0} \FloatTok{0.5}\NormalTok{) }\CommentTok{{-}{-} initial velocities}
\end{Highlighting}
\end{Shaded}
A \texttt{Config\ n} represents the state of the system, represented in
configuration-space. But, remember, Hamiltonian dynamics is about simulating the
path of the particle through \emph{phase space}. So we can convert our
configuration-space state into a phase-space state:
\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{{-}{-} source: https://github.com/mstksg/inCode/tree/master/code{-}samples/hamilton/DoublePendulum.hs\#L41{-}L42}
\OtherTok{phase0 ::} \DataTypeTok{Phase} \DecValTok{2}
\NormalTok{phase0 }\OtherTok{=}\NormalTok{ toPhase doublePendulum config0}
\end{Highlighting}
\end{Shaded}
And now we can ask for the state of our system at any amount of points in time:
\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{{-}{-} source: https://github.com/mstksg/inCode/tree/master/code{-}samples/hamilton/DoublePendulum.hs\#L44{-}L45}
\OtherTok{evolution ::}\NormalTok{ [}\DataTypeTok{Phase} \DecValTok{2}\NormalTok{]}
\NormalTok{evolution }\OtherTok{=}\NormalTok{ evolveHam\textquotesingle{} doublePendulum phase0 [}\DecValTok{0}\NormalTok{,}\FloatTok{0.1} \OperatorTok{..} \DecValTok{1}\NormalTok{]}
\end{Highlighting}
\end{Shaded}
The result there will be the state of the system at times 0, 0.01, 0.02, 0.03
\ldots{} etc.
Or, if you want to run the system step-by-step:
\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{{-}{-} source: https://github.com/mstksg/inCode/tree/master/code{-}samples/hamilton/DoublePendulum.hs\#L47{-}L48}
\OtherTok{evolution\textquotesingle{} ::}\NormalTok{ [}\DataTypeTok{Phase} \DecValTok{2}\NormalTok{]}
\NormalTok{evolution\textquotesingle{} }\OtherTok{=} \FunctionTok{iterate}\NormalTok{ (stepHam }\FloatTok{0.1}\NormalTok{ doublePendulum) phase0}
\end{Highlighting}
\end{Shaded}
And you can get the position of the coordinates as:
\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{{-}{-} source: https://github.com/mstksg/inCode/tree/master/code{-}samples/hamilton/DoublePendulum.hs\#L50{-}L51}
\OtherTok{positions ::}\NormalTok{ [}\DataTypeTok{R} \DecValTok{2}\NormalTok{]}
\NormalTok{positions }\OtherTok{=}\NormalTok{ phsPositions }\OperatorTok{<$>}\NormalTok{ evolution\textquotesingle{}}
\end{Highlighting}
\end{Shaded}
With \texttt{phsPositions\ ::\ Phase\ n\ -\textgreater{}\ R\ n}
And the position in the underlying cartesian space as:
\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{{-}{-} source: https://github.com/mstksg/inCode/tree/master/code{-}samples/hamilton/DoublePendulum.hs\#L53{-}L54}
\OtherTok{positions\textquotesingle{} ::}\NormalTok{ [}\DataTypeTok{R} \DecValTok{4}\NormalTok{]}
\NormalTok{positions\textquotesingle{} }\OtherTok{=}\NormalTok{ underlyingPos doublePendulum }\OperatorTok{<$>}\NormalTok{ positions}
\end{Highlighting}
\end{Shaded}
Where
\texttt{underlyingPos\ ::\ System\ m\ n\ -\textgreater{}\ Phase\ n\ -\textgreater{}\ R\ m}.
Let's ignore the underlying position for now, and print out now the full
progression of the system's positions:
\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{{-}{-} source: https://github.com/mstksg/inCode/tree/master/code{-}samples/hamilton/DoublePendulum.hs\#L56{-}L57}
\OtherTok{main ::} \DataTypeTok{IO}\NormalTok{ ()}
\NormalTok{main }\OtherTok{=}\NormalTok{ withRows (}\FunctionTok{take} \DecValTok{25}\NormalTok{ positions) (disp }\DecValTok{4}\NormalTok{)}
\end{Highlighting}
\end{Shaded}
(\texttt{withRows} is from \emph{hmatrix}, which treats a list of vectors as a
matrix with each vector as a row, and \texttt{disp\ 5} from \emph{hmatrix}
pretty-prints our matrix with 5 decimal places of precision)
\begin{verbatim}
L 25 2
1.0000 0.0000
0.9727 0.0800
0.8848 0.2345
0.7164 0.5129
0.4849 0.8725
0.2878 1.0648
0.1223 1.0801
-0.0165 0.9388
-0.1099 0.6400
-0.1161 0.1447
-0.0539 -0.4882
-0.0795 -0.9212
-0.1689 -1.1797
-0.2860 -1.2970
-0.4146 -1.2803
-0.5562 -1.1238
-0.7249 -0.8079
-0.8762 -0.4505
-0.9442 -0.2075
-0.9416 -0.0516
-0.8793 0.0312
-0.7596 0.0265
-0.5728 -0.1086
-0.3001 -0.4237
-0.0381 -0.6640
\end{verbatim}
Neat! We see that the first coordinate (\(\theta_1\)) starts at 1 like we asked,
and then begins decreasing and falling\ldots{} And then we see the second
coordinate (\(\theta_2\)) starting at 0 and then ``swinging'' to the right. The
\href{http://i.imgur.com/Vaaa2EC.gifv}{image the top of this post} is an
animation of such a system (albeit with \(m_2 = 1\)).
\hypertarget{two-body-system}{%
\subsubsection{Two-body system}\label{two-body-system}}
Here's one more situation where generalized coordinates describe things in a lot
nicer way than cartesian coordinates: the classic two-body problem.
Really, you can describe the state of a two-body system with only two
parameters: the distance between the two bodies, and their current angle of
rotation.
In this framework, Kepler tells us that for bodies in orbit, the distance will
grow smaller and larger again over time, and that the angle of rotation will
constantly increase\ldots and increase at a faster rate when the distance is
smaller (which is
\href{https://en.wikipedia.org/wiki/Kepler's_laws_of_planetary_motion\#Second_law}{Kepler's
second law}).
If we assume that the center of mass of the system is at
\(\langle 0, 0 \rangle\), then we can state these coordinates as
\[
\langle x_1, y_1 \rangle = \langle r_1 \cos (\theta), r_1 \sin (\theta) \rangle
\]
\[
\langle x_2, y_2 \rangle = \langle r_2 \cos (\theta), r_2 \sin (\theta) \rangle
\]
Where \(r_1 = \frac{m_2}{m_1 + m_2}\) and \(r_2 = - \frac{m_1}{m_1 + m_2}\)
(solving from the center of mass).\footnote{Alternatively, we could assume that
the halfway point (or even the first body) is always at
\(\langle 0, 0 \rangle\), but this doesn't give us as pretty of plots. The
center of mass is a nice reference point because newton's third law implies
that it remains stationary forever.}
Our potential energy function is Newton's famous
\href{https://en.wikipedia.org/wiki/Newton's_law_of_universal_gravitation}{law
of universal gravitation}:
\[
U(r, \theta) = - \frac{G m_1 m_2}{r}
\]
And, this should be enough to go for \emph{hamilton}.
``But wait,'' I hear you say. ``If we're doing a change-of-coordinate-system
into polar coordinates, don't we have to account for artifacts like centrifugal
acceleration from the fact that \(d \theta\) is non-uniform and depends on
\(r\)?''
Well, I'm glad you asked! And the answer is, nope. We don't have to account for
any weird interplay from non-uniform coordinate systems because \emph{hamilton}
arrives at the proper solution simply from the geometry of the generalized
coordinates. (And it does this using
\href{https://hackage.haskell.org/package/ad}{ad}, but more on that for a later
post!)
Anyway, here we go:
\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{{-}{-} source: https://github.com/mstksg/inCode/tree/master/code{-}samples/hamilton/TwoBody.hs\#L10{-}L42}
\OtherTok{twoBody ::} \DataTypeTok{System} \DecValTok{4} \DecValTok{2}
\NormalTok{twoBody }\OtherTok{=}\NormalTok{ mkSystem masses coordinates potential}
\KeywordTok{where}
\OtherTok{ masses ::} \DataTypeTok{R} \DecValTok{4}
\NormalTok{ masses }\OtherTok{=}\NormalTok{ vec4 }\DecValTok{10} \DecValTok{10} \DecValTok{1} \DecValTok{1}
\NormalTok{ coordinates}
\OtherTok{ ::} \DataTypeTok{Floating}\NormalTok{ a}
\OtherTok{=>} \DataTypeTok{V.Vector} \DecValTok{2}\NormalTok{ a}
\OtherTok{{-}>} \DataTypeTok{V.Vector} \DecValTok{4}\NormalTok{ a}
\NormalTok{ coordinates (}\DataTypeTok{V2}\NormalTok{ r θ) }\OtherTok{=} \DataTypeTok{V4}\NormalTok{ (r1 }\OperatorTok{*} \FunctionTok{cos}\NormalTok{ θ) (r1 }\OperatorTok{*} \FunctionTok{sin}\NormalTok{ θ)}
\NormalTok{ (r2 }\OperatorTok{*} \FunctionTok{cos}\NormalTok{ θ) (r2 }\OperatorTok{*} \FunctionTok{sin}\NormalTok{ θ)}
\KeywordTok{where}
\NormalTok{ r1 }\OtherTok{=}\NormalTok{ r }\OperatorTok{*} \DecValTok{1} \OperatorTok{/} \DecValTok{11}
\NormalTok{ r2 }\OtherTok{=} \OperatorTok{{-}}\NormalTok{ r }\OperatorTok{*} \DecValTok{10} \OperatorTok{/} \DecValTok{11}
\NormalTok{ potential}
\OtherTok{ ::} \DataTypeTok{Fractional}\NormalTok{ a}
\OtherTok{=>} \DataTypeTok{V.Vector} \DecValTok{2}\NormalTok{ a}
\OtherTok{{-}>}\NormalTok{ a}
\NormalTok{ potential (}\DataTypeTok{V2}\NormalTok{ r \_) }\OtherTok{=} \OperatorTok{{-}} \DecValTok{10} \OperatorTok{/}\NormalTok{ r }\CommentTok{{-}{-} G = 1}
\OtherTok{config0 ::} \DataTypeTok{Config} \DecValTok{2}
\NormalTok{config0 }\OtherTok{=} \DataTypeTok{Cfg}\NormalTok{ (vec2 }\DecValTok{2} \DecValTok{0}\NormalTok{) }\CommentTok{{-}{-} initial positions}
\NormalTok{ (vec2 }\DecValTok{0} \FloatTok{0.5}\NormalTok{) }\CommentTok{{-}{-} initial velocities}
\end{Highlighting}
\end{Shaded}
(we use \texttt{mkSystem} instead of \texttt{mkSystem\textquotesingle{}} because
we want to state the potential energy in terms of our generalized coordinates
\(r\) and \(\theta\))
Let's take a peek:
\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{{-}{-} source: https://github.com/mstksg/inCode/tree/master/code{-}samples/hamilton/TwoBody.hs\#L44{-}L60}
\OtherTok{phase0 ::} \DataTypeTok{Phase} \DecValTok{2}
\NormalTok{phase0 }\OtherTok{=}\NormalTok{ toPhase twoBody config0}
\OtherTok{evolution\textquotesingle{} ::}\NormalTok{ [}\DataTypeTok{Phase} \DecValTok{2}\NormalTok{]}
\NormalTok{evolution\textquotesingle{} }\OtherTok{=} \FunctionTok{iterate}\NormalTok{ (stepHam }\FloatTok{0.1}\NormalTok{ twoBody) phase0}
\OtherTok{positions ::}\NormalTok{ [}\DataTypeTok{R} \DecValTok{2}\NormalTok{]}
\NormalTok{positions }\OtherTok{=}\NormalTok{ phsPositions }\OperatorTok{<$>}\NormalTok{ evolution\textquotesingle{}}
\OtherTok{main ::} \DataTypeTok{IO}\NormalTok{ ()}
\NormalTok{main }\OtherTok{=}\NormalTok{ withRows (}\FunctionTok{take} \DecValTok{25}\NormalTok{ positions) (disp }\DecValTok{4}\NormalTok{)}
\end{Highlighting}
\end{Shaded}
\begin{verbatim}
L 25 2
2.0000 0.0000
1.9887 0.0502
1.9547 0.1015
1.8972 0.1554
1.8149 0.2133
1.7058 0.2777
1.5669 0.3523
1.3933 0.4435
1.1774 0.5647
0.9057 0.7503
0.5516 1.1413
0.2057 3.4946
0.6092 5.2275
0.9490 5.5664
1.2115 5.7386
1.4207 5.8542
1.5889 5.9424
1.7233 6.0152
1.8283 6.0785
1.9069 6.1358
1.9610 6.1892
1.9917 6.2403
1.9998 6.2904
1.9852 6.3407
1.9479 6.3923
\end{verbatim}
Neat! We see that \(r\) starts big and gets smaller, and then gets big again.
And it's clear that when \(r\) is smallest, \(\theta\) changes the fastest. Look
at it go!
Here's an animation of the same situation with some different masses:
\href{http://i.imgur.com/TDEHTcb.gifv}{\includegraphics{/img/entries/hamilton/two-body.gif}}
\hypertarget{just-you-wait}{%
\subsection{Just you wait}\label{just-you-wait}}
Now, this isn't all just useful for physics. You can state a lot of
animation/dynamics problems as motion along coordinates that aren't always
trivial. This project started, after all, as a way to simulate constant-velocity
motion along a bezier curve. (In that case, the single coordinate is the
non-uniform time parameter to the bezier curve.)
I've included more examples in the
\href{https://github.com/mstksg/hamilton\#example-app-runner}{example app
launcher} included in the library (which generated those animations you see
above), including:
\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
A spring hanging from a block sliding along a horizontal rail (a favorite of
many physics students, of course)
\item
A ball bouncing around a room, showing that you can represent bouncy walls as
potential energy functions
\item
The bezier curve example.
\end{enumerate}
Let me know in the comments if you think of any interesting systems to apply
this to, or if you have any interesting applications in physical or non-physical
ways! I'd love to hear :D
And if you're interested in the implementation using some of those Haskell
tricks I mentioned above, stay tuned :)
\hypertarget{signoff}{%
\section{Signoff}\label{signoff}}
Hi, thanks for reading! You can reach me via email at
\href{mailto:justin@jle.im}{\nolinkurl{justin@jle.im}}, or at twitter at
\href{https://twitter.com/mstk}{@mstk}! This post and all others are published
under the \href{https://creativecommons.org/licenses/by-nc-nd/3.0/}{CC-BY-NC-ND
3.0} license. Corrections and edits via pull request are welcome and encouraged
at \href{https://github.com/mstksg/inCode}{the source repository}.
If you feel inclined, or this post was particularly helpful for you, why not
consider \href{https://www.patreon.com/justinle/overview}{supporting me on
Patreon}, or a \href{bitcoin:3D7rmAYgbDnp4gp4rf22THsGt74fNucPDU}{BTC donation}?
:)
\end{document}