\documentclass[11pt,a4paper,fleqn]{scrartcl}
\usepackage{a4wide}

\usepackage[english]{babel}
\usepackage{times}

\usepackage{mathbbol}
%\usepackage[normalem]{ulem} % for \sout, \xout,...
%\usepackage{cancel}
\usepackage{amsthm}
\usepackage{amsmath}
\usepackage{mathtools}
%\usepackage{MnSymbol} % for \restriction % makes everything look wired!!!

\usepackage{paralist}
\usepackage{enumitem}
\usepackage{color}

%%% ====== Thm-Umgebungen =====
%\newtheorem{theorem}{Theorem}
%\newtheorem{proposition}{Proposition}
%\newtheorem{implication}{Implication}
%\newtheorem{lemma}[theorem]{Lemma}
%\newtheorem{corollary}[theorem]{Corollary}
%\newtheorem{definition}[theorem]{Definition}
%\newtheorem{problem}[theorem]{Problem}
%\newtheorem{solution}[theorem]{Solution}
%\newtheorem{observation}[theorem]{Observation}

%%% ====== Makros ==============

%%% ---- colors -----
\newcommand{\je}[1]{{\color{green} #1}}

\definecolor{henrik}{rgb}{.8,.3,0}
\newcommand{\cg}[1]{{\color{henrik} #1}}

\definecolor{martin}{rgb}{0,.4,1}
\newcommand{\mk}[1]{{\color{martin}#1}}

\newcommand{\rot}[1]{{\color{red}#1}}



%%% ====== math =======
\newcommand{\mc}[1]{\mathcal{#1}}
\newcommand{\mb}[1]{\mathbb{#1}}

%% ---- mathrm----
\renewcommand{\Re}{\mathrm{Re}}
\renewcommand{\Im}{\mathrm{Im}}
\newcommand{\e}{\mathrm{e}}
\renewcommand{\i}{\mathrm{i}}
\newcommand{\rmd}{\, \mathrm{d}}
\newcommand{\tr}{\mathrm{Tr}} %old
\newcommand{\Tr}{\mathrm{Tr}} %new
\newcommand{\supp}{\mathrm{supp}}
\newcommand{\dom}{\mathrm{dom}}
\newcommand{\id}{\mathrm{id}}
\newcommand{\poly}{\mathrm{poly}}
\newcommand{\landauO}{\mathrm{O}}
\newcommand{\dist}{\mathrm{dist}}
\newcommand{\spec}{\mathrm{spec}}
\newcommand{\diag}{\mathrm{diag}}
\newcommand{\GL}{\mathrm{GL}}
\newcommand{\Hom}{\mathrm{Hom}}
\newcommand{\argmin}{\mathrm{argmin}}
\newcommand{\hc}{\mathrm{h.c.}}
\renewcommand{\mod}{\mathrm{\, mod\, }}
\newcommand{\vect}{\mathrm{vec}}


%% ---- cal ----
%\newcommand{\mcL}{\mc{L}}
%\newcommand{\mcD}{\mc{D}}
%\newcommand{\mcH}{\mc{H}}
%\newcommand{\mcB}{\mc{B}}
%\newcommand{\mcA}{\mc{A}}
%\newcommand{\mcS}{\mc{S}}
%\newcommand{\mcC}{\mc{C}}
%\newcommand{\mcT}{\mc{T}}
%\newcommand{\mcK}{\mc{K}}

%% ---- mathbb -----
\newcommand{\N}{\mb{N}}
\newcommand{\Z}{\mb{Z}}
\newcommand{\R}{\mb{R}}
\newcommand{\C}{\mb{C}}
\newcommand{\1}{\mb{1}}
\renewcommand{\P}{\mb{P}}
%\newcommand{\mbL}{\mb{L}}
%\newcommand{\U}{\mb{U}}


%% ---- norms ----
\newcommand{\norm}[1]{\left\Vert #1 \right\Vert}
%\newcommand{\Norm}[1]{\left\Vert #1 \right\Vert_{1\rightarrow 1}}
%\newcommand{\iNorm}[1]{\left\Vert #1 \right\Vert_{\infty}}
%\newcommand{\iiNorm}[1]{\left\Vert #1 \right\Vert_{\infty\rightarrow\infty}}
%\newcommand{\nNorm}[1]{\lVert #1\rVert_{1\rightarrow 1}}
%\newcommand{\BNorm}[1]{\Bigl\Vert #1 \Bigr\Vert_{1\rightarrow 1}}
%\newcommand{\bnorm}[1]{\bigl\Vert #1 \bigr\Vert}
%\newcommand{\nnorm}[1]{\lVert #1\rVert}

%%% ---- Kets -----
\newcommand{\ket}[1]{\left.\left|{#1}\right.\right\rangle}
\newcommand{\ketn}[1]{| #1 \rangle}
\newcommand{\bra}[1]{\left.\left\langle{#1}\right.\right|}
\newcommand{\bran}[1]{\langle #1 |}
\newcommand{\braket}[2]{\left\langle #1 \middle| #2 \right\rangle}
\newcommand{\braketn}[2]{\langle #1 \middle| #2 \rangle}
\newcommand{\ketbra}[2]{\ket{#1} \!\! \bra{#2}}
\newcommand{\ketibra}[3]{\ket{#1}_{#2} \!\! \bra{#3}}
\newcommand{\ketbran}[2]{\ketn{#1} \! \bran{#2}}
\newcommand{\sandwich}[3]
  {\left\langle  #1 \right| #2 \left| #3 \right\rangle}

%% ---- other stuff -----
\newcommand{\kw}[1]{\frac{1}{#1}}
\newcommand{\tkw}[1]{\tfrac{1}{#1}}
\newcommand{\ad}{^\dagger}
\newcommand{\del}{\partial}
\newcommand{\ddt}{\frac{\rmd}{\rmd t}}
\newcommand{\argdot}{ \cdot }
\newcommand{\thmnl}{ $ \ $ \\ \vspace{-1\baselineskip}}
\newcommand{\ol}[1]{\overline{#1}}
\newcommand{\restr}{\restriction}
\newcommand{\bonus}[1]{\emph{Bonus question:} #1}
\renewcommand{\d}{\rm d}
\newcommand{\hint}[1]{\emph{Hint:} #1}



%%%=====================================================================
%%%========================== begin document =============================
%%%=====================================================================
\begin{document}

\begin{center}
$\ $
\vspace{-4\baselineskip}
\\
\textbf{\Large Advanced statistical mechanics II}\\
{\large (SS 16, FU Berlin)}\\ [0.5em]
\textbf{Bonus problem $0+ 1i$}\\
Due date: June 21st, 2016 
\end{center}

\subsection*{Bonus Sheet 0+1i}
\begin{enumerate}[leftmargin=-0ex]
\setcounter{enumi}{0}
%%%=====================================================================
%%%========================== 1111111111111 =============================
%%%=====================================================================
\item {\bf Evolution of phase space (1*+1*+1*+1*+1*+10*+5*+0*=20*)}\\
How does phase space evolve in time? 
There is no better way than to see for yourself in a cute simulation!
In this exercise we will study the evolution of phase space $\mathcal P=\{(q,p)\}=\mathbb R^{\times 2}$ of a single particle confined to a one dimensional potential.
\begin{enumerate}
\item \label{eq:dS}
 Write down the (classical) Hamiltonian  and the equations of motion of a particle subject to a potential $V(q)$.
\end{enumerate}
From now on, we set the mass of the particle to $1$.
\begin{enumerate}[resume]
\item Consider $V(q)=\frac{1}{2}q^2(\frac{1}{2}q^2-1)$. Find the function $\phi((q,p),t)\in \mathbb R^{2\times 1}\rightarrow\mathbb R^{2}$ such that 
\begin{align}
\begin{bmatrix} 
\dot{p}\\
\dot{q}  
\end{bmatrix}
=\phi((q,p),t)
\end{align}
What does $\phi$ actually depend on? What is needed to know the trajectory of the particle?
\end{enumerate}
We will now try to explain how software solving differential equations works. An ordinary differential equations (ODE) is a problems of the following form 
  \begin{align}
 \text{\noindent "Find a function }y: I\rightarrow \mathbb R\text{ s.t. }\forall t\in I\quad F(y,dy/dt,\ldots,d^ny/dt^n,t)=0"\;.
  \end{align}
 For example, $F$ could be $d^2y/dt^2+\omega^2 y=0$.
 It turns out that another form is more handy and in numerics one considers the form
\begin{align}
  \frac{dx}{dt}=f(x,t)\quad x\in \mathbb R^N \label{eq:ode_form}
\end{align}
Let us motivate why the formulation with a \emph{single} derivative is not less general.
\begin{enumerate}[resume]
\item
 As an example consider the Van der Pol oscillator with time-dependent driving $\theta(t)=3cos(t)$
 \begin{align}
   \frac{d^2x}{dt^2}-\mu(1-x^2) \frac{dx}{dt}+x+\theta(t)=0\;.\label{eq:vdP}
 \end{align}
Introduce redundant variables $p=\frac{dx}{dt}, q=x, \lambda =t$ and find the vector notation  in $y=(p,q,\lambda)$ the equivalent first-order formulation of \eqref{eq:vdP} reading $\dot y=f(y)$ for some $f$ that you should write out explicitly. 
Argue that if you know a vector-valued function $g=(g_1,g_2,g_3)$ satisfying $\dot g=f(g)$ in a given region then its $g_2$, which is the entry corresponding to the coordinate $q$, solves \eqref{eq:vdP}.

\hint{$d\lambda/dt=1$}

\item Use the fact that we may consider \eqref{eq:ode_form} to argue that an ODE of order $n$ needs an initial condition consisting of $n$ values.
Write down the ODE of a harmonic oscillator and show how this applies to this example, i.e. how many constants do you need for the solution to be uniquely specified? 
How many initial values does an ideal gas of $N_A$-many classical particles need?
\end{enumerate}
To solve an ODE we discretise the time and use Taylor's expansion\footnote{You may want to read up on solvers constructed from Trotterization - those are useful in hybrid Monte-Carlo methods for molecular dynamics.}. The simplest one is the Euler rule.
\begin{enumerate}[resume]
\item
Consider a function which for $N\in \mathbb Z_>$ takes an initial value $r_0=(p_0,q_0)$ and time step $\Delta t=1/N$ as input and gives back values $\{r_0,r_1,\ldots r_N\}$.
Assuming that for all $\Delta t$ we have
\begin{align}
  r_k=r_{k-1}+f(r_{k-1},(k-1)\Delta t)\Delta t\quad k\in[N]
\end{align}
argue that $r$ satisfies $\dot r(t)=f(r(t),t)$ for $t\in[0,1]$ in the limit $\Delta t\rightarrow 0$.
\end{enumerate}


With the exercise comes the Octave script BonusOctave.m which can be easily used on the computers at the university. 
It integrates the equations of motion of a particle in a Mexican-hat potential.
In order to model the flow of the phase space in time we discretise a rectangle of size $[-0.1,0.1]^{\times 2}$ into a grid and look at the evolution of of phase space points on the nodes. 
By integrating the equations of motion for each initial condition we can see how the phase space evolves as in the gif on the tutorial page for this exercise.

In the next exercise you will be asked to produce your own output by changing the potential.  
To use the script 'BonusOctave.m' provided with this exercise do the following.
 Create a directory for this numerical project, say 'Bonus' and 'Bonus/output' where you we will keep the snapshots of the simulation.
 Move the script 'BonusOctave' to the directory 'Bonus'.
 Use the console to cd into the directory 'Bonus' and type 'octave' in the console to run Octave. In octave run 'BonusOctave()' to start the simulation.
The integration should take no longer than a minute, however rendering all the snapshots may take a bit. You should see 100 pictures enumerated from 00001.png to 00100.png.

To create a gif run: "convert -delay 10 -loop 0 output/*.png ThisIsFun.gif"
\begin{enumerate}[resume]
\item Design your own potential $V(q)$, implement it appropriately as the function 'f(x)=\dot{x}' at the end of 'BonusOctave.m' and explore what is possible with this simple script!

Create a pdf with a description of what you did, what worked and what did not.
Put in representable snapshots that illustrate your observations.

Possible questions to discuss:
Do you think it is a coincidence that the box is being sheared in the initial stages quite deterministically into a rhombus?\\
For how long can you tell what was the initial phase space configuration?\\
How is the phase space populated restricted to an energy window?\\
Is there a difference, whether you evolve forward or backwards in time?\\
Can you observe a recurrence to the initial configuration?\\
Can you interpret the phase space evolution in terms of an ideal gas?\\
What would you need to model 30 particles in a box with van der Waals interactions?

\item Modify the script 'BonusOctave.m' so that in your animation you see the particles in real space.
For example, on each frame show two plots: on the left an ensemble of beads travelling on a curved wire (i.e. plot the position of each initial conditions and its potential energy) and on the right the phase space evolution.

\item Can we use your results to present it in class?
\end{enumerate}
\subsection*{Regulations for the bonus problem}
\begin{compactitem}
\item In order to redistribute your workload you have one month to finish this bonus exercise
\item It's a good occasion to get lots of points for your own satisfaction! 
\item Preferably use latex to show what you did and gifs for the in-class presentation 
\item You can work with the people you like, but please provide at least one animation for each person. 
\item Last, but not least: Enjoy! :)
\end{compactitem}


\end{enumerate}
\end{document}
