\documentclass[11pt]{article}
\usepackage{fullpage,amsmath,amsfonts,mathpazo,microtype,nicefrac,algorithm2e,graphicx}

% Set-up for hypertext references
\usepackage{hyperref,color,textcomp}
\definecolor{webgreen}{rgb}{0,.35,0}
\definecolor{webbrown}{rgb}{.6,0,0}
\definecolor{RoyalBlue}{rgb}{0,0,0.9}
\hypersetup{
   colorlinks=true, linktocpage=true, pdfstartpage=3, pdfstartview=FitV,
   breaklinks=true, pdfpagemode=UseNone, pageanchor=true, pdfpagemode=UseOutlines,
   plainpages=false, bookmarksnumbered, bookmarksopen=true, bookmarksopenlevel=1,
   hypertexnames=true, pdfhighlight=/O,
   urlcolor=webbrown, linkcolor=RoyalBlue, citecolor=webgreen,
   pdfauthor={Chris H. Rycroft},
   pdfsubject={Harvard AM225 (Spring 2021)},
   pdfkeywords={},
   pdfcreator={pdfLaTeX},
   pdfproducer={LaTeX with hyperref}
}
\hypersetup{pdftitle={AM225: Assignment 1}}

\SetKwProg{Fn}{}{}{}\SetKwFunction{FRecurs}{void SudokuSolve}

% Macro definitions
\newcommand{\N}{\mathbb{N}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\p}{\partial}
\renewcommand{\vec}[1]{\mathbf{#1}}
\newcommand{\vx}{\vec{x}}
\newcommand{\vv}{\vec{v}}
\newcommand{\vp}{\vec{p}}
\newcommand{\rs}{\sqrt{6}}
\newcommand{\Atol}{\textit{Atol}}
\newcommand{\Rtol}{\textit{Rtol}}
\newcommand{\Trans}{\mathsf{T}}

\DeclareMathOperator{\signop}{sign}

\begin{document}
\section*{AM225: Assignment 2 (due 5~PM, March 4)}
\subsection*{Part I: ODE solution methods}
\noindent \textit{Complete at least \textbf{one out of two} problems in this
section. If you submit answers for both, your grade will be calculated using
the best score. Note: question 2 is harder than question 1.}
\begin{enumerate}
  \item \textbf{Adaptive integration with a First Same As Last (FSAL) scheme.}
      Write an adaptive Runge--Kutta integration scheme to solve an arbitrary
      ODE problem $y'=f(x,y)$ for $y(x)\in \R^n$, using the five step FSAL
      scheme described in the lectures. Your code should have the following
      properties:
      \begin{itemize}
	\item Use the fourth-order solution $y_1$ to step forward, and use the
	  third-order $\hat{y}_1$ for step size selection.
	\item Use the step size selection procedure described in the slides.
	  The parameters $\textit{fac}$, $\textit{facmax}$, and
	  $\textit{facmin}$ from the slide can be used. For the tolerance, you
	  can assume $\Atol$ and $\Rtol$ are the same for all
	  components, rather than being specified on a per-component basis.
	  Use an initial step size of $h=0.01$.
	\item Make your program count the number of function evaluations.
	  Your program should be parsimonious in the number of function
	  evaluations it performs, by reusing $k_5$ from a succesful step to
	  be $k_1$ of the next step, and by retaining $k_1$ when a step is
	  rejected.
      \end{itemize}
      Once your method is working, complete the following two test problems.
    \begin{enumerate}
      \item Test your program on the Brusselator test problem using
	$\Atol=\Rtol=\lambda$. By trying a range of $\lambda$ from $10^{-3}$ to
	$10^{-13}$ make a precision--work plot as in the lectures. In the
	program files, you will find the corresponding precision--work data for
	the four low-order methods considered in the lectures, and you may
	overlay these results on your plot to compare them.
      \item Extend your code so that it does third-order dense output with
	Hermite interpolation at regular intervals. Test your code using the
	two-component system
	\begin{equation}
	  y_1'=-xy_2, \qquad y_2' = xy_1
	\end{equation}
	with initial conditions $y_1(0)=1, y_2(0)=0$. This problem has the
	exact solution \smash{$y_1^\text{exact}(x)=\cos \frac{x^2}{2},
	y_2^\text{exact}=\sin \frac{x^2}{2}$}. Simulate to $x=8$ using
	$\lambda=3\times 10^{-3}$, saving dense output at intervals of
	\smash{$\frac{8}{1200}$}. Plot the numerically computed solutions
	$y_1^\text{num}$ and $y_2^\text{num}$, showing the integration steps as
	points, and the dense output as lines. Make a second plot showing
	$y_1^\text{num}-y_1^\text{exact}$ and
	$y_2^\text{num}-y_2^\text{exact}$.
      \item \textbf{Optional.} The method is not sensitive to the initial step
	size choice of $h=0.01$, since the adaptive procedure will
	automatically adjust it. However, Hairer \textit{et al.}~describe
	an algorithm for estimating the initial timestep. Extend your
	method to implement this.
    \end{enumerate}
  \item \textbf{A high-order adaptive integrator using Richardson
    extrapolation.} Write an adaptive Runge--Kutta integration scheme by
    applying Richardson extrapolation to the fifth-order Cash--Karp
    scheme\footnote{For the purposes of this question, you can ignore the lower
    order Cash--Karp formulae, since here the aim is to use Richardson
    extrapolation for step size selection.} in the lectures, thereby obtaining
    a sixth-order method. Starting from $y_0$, let $y_1$ and $y_2$ be
    Cash--Karp steps with size \smash{$\frac{h}{2}$}, and $w$ be a
    Cash--Karp step of size $h$. Define the sixth-order solutions
    \begin{equation}
      \hat{y}_1 = y_1 + \frac{y_2-w}{(2^p-1)2}, \qquad \hat{y}_2 = y_2 + \frac{y_2 -w}{2^p-1}.
    \end{equation}
    Your program should use the same step size selection procedure as from
    Question 1.\footnote{Since Richardson extrapolation requires taking two
    timesteps of size $h/2$, you may encounter NaNs for a large choice of $h$.
    You code should reject that step and try again with $h\times \textit{facmin}$.} It
    should count the number of function evaluations and be as parsimonious as
    possible. Use $y_2-\hat{y}_2$ for step size selection, and use $\hat{y}_2$
    to advance forward in $x$.
    \begin{enumerate}
      \item Repeat Question 1(a) for this method.
      \item Extend your code so that it computes dense output at regular
	intervals, based on quintic polynomial interpolation using $y_0$,
	$f(x_0,y_0)$, $\hat{y}_1$, $f(x_0+h,y_1)$, $\hat{y}_2$, and
	$f(x_0+2h,\hat{y}_2)$.\footnote{Note that the derivative $f(x_0+h,y_1)$
	can be used. It is not necessary to evaluate $f(x_0+h,\hat{y}_1)$.}
	Repeat the two-component test from Question 1(b).
    \end{enumerate}
\end{enumerate}

\subsection*{Part II: ODE applications and analysis}
\noindent \textit{Complete at least \textbf{three out of five} problems in this
section. If you submit answers for more, your grade will be calculated using
the three best scores.}

\begin{enumerate}
  \setcounter{enumi}{2}
  \item \textbf{Order condition trees.} Write a program to enumerate all trees
    of a given order. Provide a list of the number of trees up to order
    15.\footnote{In the lecture slides you will find the number of trees up to
    order 10, which you can use to check your solutions.} Extend your program
    so that it can visualize the trees in some format of your choice, and use
    it to show all trees of order 7.
  \item \textbf{Error analysis of a Richardson extrapolation scheme.}
    \begin{enumerate}
      \item Show that Richardson extrapolation applied to the second-order
	Ralston method can be reformulated as a five-step, third-order
	Runge--Kutta method, and find its Butcher tableau.
      \item The third-order Heun method has Butcher tableau
	\[
	\begin{array}{c|ccc}
	  0 & & & \\
	  \nicefrac{1}{3} & \nicefrac{1}{3} & & \\
	  \nicefrac{2}{3} & 0 & \nicefrac{2}{3} & \\
	  \hline
	  & \nicefrac{1}{4} & 0 & \nicefrac{3}{4} \\
	\end{array}
        \]
	For both the Heun method, and your method from part (a), determine the
	error coefficients $e(t)$ for all trees $t$ of order 4, reporting your
	answers as rational numbers.
      \item Show that one of the methods has universally smaller error
	magnitudes $|e(t)|$ than the other. Once the difference in the number
	of function evaluations is taken into account, will that method be
	better for practical calculations?
    \end{enumerate}
  \item \textbf{A generalized Kuramoto model.}\footnote{This question was
    suggested by Nick Boffi (\texttt{boffi@g.harvard.edu}), and could be the
    basis for a final project.} A recent paper by O'Keeffe
    \textit{et al.}\footnote{K.~P.~O'Keeffe, H.~Hong, and S.~H.~Strogatz,
    \textit{Oscillators that sync and swarm}, Nat.\ Commun.\ \textbf{8}, 1504
    (2017).
    \href{http://dx.doi.org/10.1038/s41467-017-01190-3}{\texttt{doi:10.1038/s41467-017-01190-3}}}
    explores a model for swarming and synchronization behavior. In the model,
    we consider $N$ agents with positions $\vx_i(t)$ and internal phases
    $\theta_i(t)$, which move according to the differential equations
    \begin{align}
      \dot{\vx}_i &= \vv_i + \frac{1}{N} \left[ \sum_{j\ne i}^N \frac{\vx_j-\vx_i}{|\vx_j-\vx_i|} (A+J\cos (\theta_j-\theta_i))
      - B\frac{\vx_j-\vx_i}{|\vx_j-\vx_i|^2} \right], \label{eq:swarm1} \\
      \dot{\theta}_i &= \omega_i + \frac{K}{N} \sum_{j\ne i}^N \frac{\sin(\theta_j -\theta_i)}{|\vx_j-\vx_i|} \label{eq:swarm2}
    \end{align}
    where $J$ and $K$ are constants, and $\vv_i$ and $\omega_i$ can be
    individually controlled for each agent. By rescaling time and space, we can
    restrict attention to the case when $A=B=1$.
    \begin{enumerate}
      \item By making use of your favorite adaptive integrator with dense
	output\footnote{You could use your code from Q1 or Q2. You could use
	the DOP853 implementation found in the course example codes.} solve
	Eqs.~\ref{eq:swarm1} \& \ref{eq:swarm2} using $N=1250$ agents. Set
	$\vv_i=\omega_i=0$ for all agents. Simulate from $t=0$ to $t=200$, and
	use dense output to save the positions at $n$ equally-spaced intervals,
	where $N\ge401$. Use $\textit{Atol}=\textit{Rtol}=10^{-6}$ in your
	adaptive integration routine.

	Use initial conditions of random positions in the unit disk,
	$\|\vx\|\le 1$, and random phases over $[0,2\pi)$. Visualize the agents
	as dots that are colored according to their phase. A suggested color
	palette is
	\[
	(R,G,B)=(f(\theta),f(\theta-\nicefrac{2\pi}{3}),
	f(\theta+\nicefrac{2\pi}{3}))
	\]
	where $f(\theta)=0.45(1+\cos\theta)$. Simulate the model with the
	following parameters:
	\begin{enumerate}
	  \item $J=0.5,K=0.5$,
	  \item $J=0.3,K=-0.2$,
	  \item $J=1,K=-0.2$.
	\end{enumerate}
	For each case, state the total number of timesteps taken. Either
	\begin{itemize}
	  \item include snapshots after $t=10,20,50,200$,
	  \item or make a movie of the snapshots.
	\end{itemize}
      \item Simulate at least one possible variation. Examples include: (i)
	changing $\vv_i$ and $\omega_i$, (ii) simulating two systems to steady
	state and then making a new initial condition with both superimposed,
	and (iii) implementing the method in 3D.\footnote{See O'Keeffe
	\textit{et al.}~to see how they alter the strengths of the terms
	in 3D.}
      \item \textbf{Optional.} Extend your code to calculate right hand sides
	of Eqs.~\ref{eq:swarm1} \& \ref{eq:swarm2} using OpenMP. Note that the
	influence of actor A on actor B is equal and opposite to the influence
	of actor B on actor A. Structure your code so that it only considers
	each pair once.
    \end{enumerate}
  \item \textbf{Symplectic integration for galactic dynamics.} The
    following fifth-order IRK method due to Geng is symplectic, meaning
    that it exactly preserves the Hamiltonian $H(p,q)$ for a Hamiltonian
    system:
    \def\arraystretch{2.5}
    \[
    \begin{array}{c|ccc}
      \dfrac{4-\sqrt{6}}{10} & \dfrac{16-\rs}{72} & \dfrac{328-167\rs}{1800} & \dfrac{-2+3\rs}{450} \\
      \dfrac{4+\sqrt{6}}{10} & \dfrac{328+167\rs}{1800} & \dfrac{16+\rs}{72} & \dfrac{-2 -3\rs}{450}\\
      1 & \dfrac{85-10\rs}{180} & \dfrac{85+10\rs}{180} & \dfrac{1}{18} \\
      \hline
      & \dfrac{16-\rs}{36} & \dfrac{16+\rs}{36} & \dfrac{1}{9} \\
    \end{array}
    \]
    \def\arraystretch{1}
    A simple model for the movement of star in a galaxy is described by the
    Hamiltonian
    \begin{equation}
      H(\vec{p},\vec{q})= \frac{p_1^2 + p_2^2 + p_3^2}{2} + \Omega (p_1q_2 -p_2q_1) + V(\vec{q}),
    \end{equation}
    where the star's position is $\vec{q}(t)=(q_1(t),q_2(t),q_3(t))$ and
    its momentum is $\vec{p}(t)=(p_1(t),$ $p_2(t), p_3(t))$. Here, $\Omega$ is the
    galaxy's velocity, and $V$ is the gravitational potential, which is
    approximated as
    \begin{equation}
      V(\vec{q}) = A \log\left(C + \frac{q_1^2}{a^2} + \frac{q_2^2}{b^2} +
      \frac{q_3^2}{c^2}\right).
    \end{equation}
    We use non-dimensionalized parameters $a=1.25$, $b=1$, $c=0.75$, $A=1$,
    $C=1$, $\Omega=0.25$. The Hamiltonian differential equation system is
    given by
    \begin{equation}
      \dot{p}_i = -\frac{\p H}{\p q_i}, \qquad \dot{q}_i = \frac{\p H}{\p p_i}
    \end{equation}
    for $i=1,2,3$. Initial conditions are given by $q_2(0)=q_3(0)=p_1(0)=0$,
    $q_1(0)=2.5$, and $p_3(0)=0.2$. The remaining momentum coordinate is
    chosen to be the larger of the two roots that yields $H=2$.
    \begin{enumerate}
      \item Implement Geng's method, and test it on an ODE of your choice
	to verify that it is fifth-order accurate.\footnote{You do not need
	to test the method on a symplectic ODE system.} Make a
	convergence plot demonstrating fifth-order accuracy.
      \item Simulate the galaxy ODE system up to $t=2000$ using a step size
	of $\nicefrac{1}{20}$ and make a 3D plot of the trajectory. Plot
	the Hamiltonian up to $t=2000$.
      \item Simulate up $t=10^5$ and make a
	\href{https://en.wikipedia.org/wiki/Poincare_map}{Poincar\'e map}
	by tracking all intersections with the half-plane $q_1>0$, $q_2=0$,
	where $\dot{q}_2>0$. You can find the intersection points by
	approximating the trajectory as a linear segments between
	successive timesteps.
    \end{enumerate}
  \item \textbf{Integrating ODEs with discontinuities.} Consider the
    two-component ODE system for functions $x(t)$ and $y(t)$ given by
    \begin{equation}
      \frac{dx}{dt}=
      \begin{cases}
	0 & \qquad \text{if $|x| \ge |y|$,} \\
	-y & \qquad \text{if $|x| < |y|$,}
      \end{cases}
      \label{eq:dx_dis}
    \end{equation}
    and
    \begin{equation}
      \frac{dy}{dt}=
      \begin{cases}
	x & \qquad \text{if $|x| \ge |y|$,} \\
	0 & \qquad \text{if $|x| < |y|$.}
      \end{cases}
      \label{eq:dy_dis}
    \end{equation}
    Use the initial condition $x(0)=1$ and $y(0)=0$.
    \begin{enumerate}
      \item Calculate the analytical solutions of $x(t)$ and $y(t)$. Show
	that they are periodic, and find the period.
      \item Simulate the ODE system in Eqs.~\ref{eq:dx_dis} and
	\ref{eq:dy_dis} to $t=48+e^{-1}$, using the classic fixed-step
	fourth-order Runge--Kutta (RK4) method. Make a work--precision plot
	using a range of total step numbers from $10^3$ to $10^7$. Your
	calculation of precision should be based on the difference between the
	numerical solution and the exact solution from part (a). Is the
	convergence data consistent with RK4 being fourth-order accurate? If
	not, why not?
      \item Repeat part (b) with your favorite adaptive integrator, using
	$\Rtol=0$, and a range of absolute tolerances of $\Atol\in
	[10^{-12},10^{-2}]$. Overlay the results on the work--precision
	plot from part (b).
      \item \textbf{Optional.} Consider the variant ODE system
	\begin{equation}
	  \frac{dx}{dt}=
	  \begin{cases}
	    0 & \qquad \text{if $|x| \ge |y|$,} \\
	    -\signop(y) & \qquad \text{if $|x| < |y|$,}
	  \end{cases}
	  \label{eq:dx_dis2}
	\end{equation}
	and
	\begin{equation}
	  \frac{dy}{dt}=
	  \begin{cases}
	    \signop(x) & \qquad \text{if $|x| \ge |y|$,} \\
	    0 & \qquad \text{if $|x| < |y|,$}
	  \end{cases}
	  \label{eq:dy_dis2}
	\end{equation}
	with initial conditions $x(0)=1$ and $y(0)=0$. Here
	\begin{equation}
	  \signop(x)=
	  \begin{cases}
	    1 & \qquad \text{if $x>0$,} \\
	    0 & \qquad \text{if $x=0$,} \\
	    -1 & \qquad \text{if $x<0$.}
	  \end{cases}
	\end{equation}
	Find the exact solution for this ODE system. Repeat the convergence
	analysis from parts (b) and (c), and compare the work--precision plots.
    \end{enumerate}
\end{enumerate}
\end{document}
