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

% 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 4}}
\usepackage[font=footnotesize,skip=3pt]{caption}

\graphicspath{ {am225_hw4_files/} {./} }

% 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{\calC}{\mathcal{C}}
\newcommand{\calF}{\mathcal{F}}
\newcommand{\Atol}{\textit{Atol}}
\newcommand{\Rtol}{\textit{Rtol}}
\newcommand{\Trans}{\mathsf{T}}
\newcommand{\sep}{\,:\,}

\DeclareMathOperator{\atantwo}{atan2}

\renewcommand{\topfraction}{0.9}
\renewcommand{\bottomfraction}{0.8}
\renewcommand{\textfraction}{0.2}
\renewcommand{\floatpagefraction}{0.7}
\renewcommand{\dblfloatpagefraction}{0.7}

\begin{document}
\section*{AM225: Assignment 4 (due 5~PM, April 16)}
\noindent \textit{Complete \textbf{at least two} questions. If you submit
answers for more, your grade will be calculated using the best scores.}

\begin{enumerate}
  \item \textbf{A previous question.} Complete a question from HW3 that you
    have not already done. In particular, question 3 from HW3 on the
    finite-element method is recommended.
  \item \textbf{Fitting a square peg in a round hole.} Consider two coordinate
    systems $\vv=(v,w)$ and $\vx=(x,y)$. In the $\vv$ system, introduce the
    circle of unit radius, $\Omega=\{ \vv\in\R \sep \|\vv\|_2<1\}$. In the
    $\vx$ system, introduce the square $S=(-1,1)^2$. The mapping $\vv=f(\vx)$
    defined by
    \begin{equation}
      v=x\sqrt{1-\frac{y^2}{2}}, \qquad w=y\sqrt{1-\frac{x^2}{2}}
    \end{equation}
    is a differentiable map from $S$ to $\Omega$ (Fig.~\ref{fig:squircle}).
    Consider the Poisson problem
    \begin{equation}
      - \nabla^2 u = f
    \end{equation}
    on the circle $\Omega$ with Dirichlet condition $u(\vv)=0$ for $\vv\in \p
    \Omega$. Introduce an $N \times N$ grid of squares on $S$, each with side
    length $h=2/N$. Let $\phi_i$ be a set of finite element basis functions on
    $S$. You can use the bilinear elements corresponding to the nodal basis at
    square corners.\footnote{You can use an alternative basis if you prefer,
    such as the $C^1$ Bogner--Fox--Schmit element.} Since the elements on the
    boundary can be neglected, this gives $(N-1)^2$ basis functions in total.
    On the circle, define basis functions $\psi_i$ via the mapping from $S$,
    such that
    \begin{equation}
      \psi_i(\vv) = \phi_i(f^{-1}(\vv)).
    \end{equation}
    Consider representing your solution as
    \begin{equation}
      u(v,w)= \sum_i u_i \psi(v,w)
    \end{equation}
    and solve the PDE problem using the Ritz--Galerkin method. Use a source
    term of $f(v,w)=e^{-v}(3+(v-4)v+w^2)$, which gives the analytical solution
    $u(v,w)=(1-v^2-w^2)e^{-v}$. For this problem, the stiffness and mass matrix
    calculations will vary from element to element, and you can compute them
    by pulling back the integrals from $\Omega$ to $S$. For example
    \begin{align}
      a(\psi_i,\psi_j) &= \int_\Omega \nabla_\vv \psi_i \cdot \nabla_\vv \psi_j \, d\vv \nonumber \\
      &= \int_S (D^{-1} \nabla_\vx \phi_i) \cdot (D^{-1} \nabla_\vx \phi_j)
      (\det D) d\vx \label{eq:fe_pull_back}
    \end{align}
    where $D = \p \vv/\p \vx$ is the Jacobian of the mapping. In general these
    integrals will not analytically solvable, but they can be performed
    accurately and efficiently using Gaussian quadrature. In the class example
    files, you will find a demo program \texttt{3a\_f\_element/q2d\_test.cc},
    which demonstrates integrating a function on a square using Gaussian
    quadrature.\footnote{Note that because each integral like
    Eq.~\ref{eq:fe_pull_back} is only done over a small patch, the integrand
    should be well-approximated by a low-order Taylor series, and thus not
    many quadrature points are require to achieve very high accuracy.}

    Use your code for a variety of choices of $N$ to calculate the $L_2$ error
    between the numerical solution and analytical solution, and determine the
    rate of convergence.
\begin{figure}
  \begin{center}
    \includegraphics[width=0.6\textwidth]{squircle_fig}
  \end{center}
  \vspace{-1em}
  \caption{The mapping $f$ from the square $S$ to the circle $\Omega$ used in question 2.\label{fig:squircle}}
\end{figure}
  \item \textbf{A generalization of the Lax--Wendroff scheme.} Consider the
    hyperbolic conservation equation
    \begin{equation}
      q_t + [A(x) q]_x =0 \label{eq:cons}
    \end{equation}
    for a function on $q(x,t)$ on the periodic interval $[0,2\pi)$. Let
    \smash{$A(x)=2+\tfrac{4}{3}\sin x$}. Following the finite volume approach,
    divide the intervals into $m$ domains $\calC_i$ of length
    $h=\tfrac{2\pi}{m}$, for $i=\{0,1,\ldots,m-1\}$. Let $Q_i^n \approx q(
    (i+\nicefrac12)h, n\Delta t)$ be the discretized solution at the center of
    each $\calC_i$. The generalized Lax--Wendroff scheme for this equation is
    given by
    \begin{equation}
      Q^{n+1}_i = Q_i - \frac{\Delta t}{h} \left[ \calF_{i+1/2}^n - \calF_{i-1/2}^n\right]
    \end{equation}
    where the fluxes are
    \begin{equation}
      \calF_{i-1/2}^n = \frac{A_{i-1} Q_{i-1}^n + A_i Q_i^n}{2} - \frac{A_{i-1/2}\Delta t}{2h} \left[A_iQ^n_i - A_{i-1}Q_{i-1}^n\right]. \label{eq:lw}
    \end{equation}
    Here, $A_i = A( (i+\nicefrac12)h)$ and $A_{i-1/2}=A(ih)$.
    \begin{enumerate}
      \item By the considering the characteristics, or otherwise, show that the
	mathematical solution to Eq.~\ref{eq:cons} is time-periodic, so that
	$q(x,t+T)=q(x,t)$ for some $T>0$. Determine $T$.\footnote{This is
	not as straightforward as it may seem. You may need to use a symbolic
	solver such as Mathematica.}
      \item The CFL condition requires that \smash{$\Delta t \le \tfrac{h}{c}$}
	for stability. What is $c$ in this case?
      \item Implement Eq.~\ref{eq:lw} and set $\Delta t = \tfrac{h}{3c}$. Use
	the initial condition
	\begin{equation}
	  q(x,0) = \exp\left(\sin x + \tfrac{1}{2}\sin 4x\right). \label{eq:lwic}
	\end{equation}
	For $m=512$, plot snapshots of the solution for
	$t=0,\tfrac{T}{4},\tfrac{T}{2},\tfrac{3T}{4},T$.\footnote{Since
	multiples of $\Delta t$ do not exactly match these snapshot times, you
	may need to make a small adjustment to the timestep.}
      \item By considering a range of $m$ (\textit{e.g.}~256 and upward) with
	the initial condition in Eq.~\ref{eq:lwic} calculate the $L_2$ norm
	between the numerical solution at $t=T$ and the exact answer. Determine
	the order of convergence.\footnote{When determining the order of
	convergence, you are interested in the asymptotic properites of error
	as $m$ gets large. You can ignore initial transients in error.}
      \item Repeat parts (c) and (d) for the initial condition
	\begin{equation}
	  q(x,0) = \max\{ \tfrac{\pi}{2}-|x-\pi|,0\}.
	\end{equation}
    \end{enumerate}
\end{enumerate}
\end{document}
