text:
\section{Numerical examples}
\label{sec:numerical_examples}

In this section we demonstrate the accuracy and robustness of the proposed numerical scheme on several standard
test cases for the Burgers equation \cref{eqn:burgers}, the shallow water equations \cref{eqn:shllw}, and
the compressible Euler equations \cref{eqn:euler}.

\subsection{\it Burgers equation}
\label{subsec:burgers_example}
We begin by applying the proposed method to the Burgers equation \cref{eqn:burgers} with periodic boundary conditions
on $x \in [0,1]$ and the smooth initial data:
\begin{equation}
\label{eqn:burgers_example}
  q(t=0,x) = q_0(x) = \sin\left(2 \pi x \right).
\end{equation}
This initial condition forms a stationary shock at $x=0.5$ at time
\begin{equation}
\label{eqn:burgers_example_shock_time}
t_{\text{shock}} =  \left[ \max_{\xi \in \reals} \left\{ -2\pi\cos\left( 2 \pi \xi \right) \right\} \right]^{-1}
 = \left(2\pi \right)^{-1} \approx 0.15915494309.
\end{equation}
In this example, the only limiter that is active is our variant of the Krivodonova \cite{article:Kriv07} correction step limiter
(see \cref{subsec:nonoscillatory}); and therefore, we can use this example to isolate its efficacy.
We compare two simulations that are identical in every way except one: (1)
in the first simulation no limiter is applied, and (2) in the second simulation the correction step limiter is used.
The results with $\melems = 100$, $\morder=4$, and at time $t=5(4\pi)^{-1} \approx 0.39788735773$ are shown
in \cref{fig:burgers}. The individual panels show (a) the solution with no limiters,
(b) a zoomed-in version of this solution, (c) the solution with the correction step limiter, and
(d) a zoomed-in version of this solution. In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution.

From these results we can draw some conclusions: (1) unsurprisingly, without a limiter the numerical solution is oscillatory,
producing unphysical overshoots and undershoots near the shock location; (2) the correction step limiter, which is applied only once per time-step (i.e., it is not applied in any of the prediction step iterations), successfully damps out the unphysical oscillations; and (3) the correction step limiter is not overly aggressive in that the damping effect is apparent only in the two elements on either side of the stationary shock.

\begin{figure}[!t]
\centering
\begin{tabular}{cc}
(a)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_burgers_no_limiter.pdf} &
(b)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_burgers_zoom_no_limiter.pdf} \\
(c)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_burgers.pdf} &
(d)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_burgers_zoom.pdf}
\end{tabular}
\caption{Numerical and exact solution of Burgers equation \cref{eqn:burgers}.
Shown are solutions computed with a smooth initial condition \cref{eqn:burgers_example} that forms a shock at
time \cref{eqn:burgers_example_shock_time}, periodic boundary conditions,
$\melems = 100$, and $\morder=4$ (i.e., fourth-order scheme with polynomial degree $\mdeg=3$). The solutions are shown at time $t=5(4\pi)^{-1} \approx 0.39788735773$, well after the shock has formed.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution. The individual panels show (a) the solution with no limiters,
(b) a zoomed-in version of this solution, (c) the solution with the correction step limiter, and
(d) a zoomed-in version of this solution.
 \label{fig:burgers}}
\end{figure}

\subsection{\it Shallow water equations}
We next consider solving the shallow water equations \cref{eqn:shllw} with various initial conditions that
demonstrate the accuracy of the scheme and the efficacy of the proposed limiting strategies.
In all examples the gravitational constant is taken as $g=1$ and and the limiter parameter is taken as
$\varepsilon = 10^{-14}$.

\subsubsection{Convergence test}
In order to verify the order of accuracy of the proposed scheme we make use of the so-called
{\it method of manufactured solutions} (e.g., see \cite{book:KnSa02}), where we prescribe a ``solution'',
and then add a source term to our original PDE that guarantees that this ``solution'' satisfies the PDE with
the additional source term. The additional source term will depend explicitly on time and space, but not on
the conserved variables $q$.

We consider an example where the manufactured solution is given by
\begin{equation}
\label{eqn:meth_man_solns}
h_{\text{MS}}(t,x) := 1 + \frac{1}{2} \sin \left(\pi (x-t) \right), \quad
u_{\text{MS}}(t,x) := \cos\left(2 \pi (x-2t) \right),
\end{equation}
on the domain $[-1,1]$ with periodic boundary conditions.
The shallow water equations, with an additional source term to guarantee that \cref{eqn:meth_man_solns}
is indeed a solution, can be written as follows:
\begin{equation}
\label{eqn:forced_shllw}
\begin{bmatrix}[1.5] h \\ hu \end{bmatrix}_{,t} +
\begin{bmatrix}[1.5] hu \\ hu^2 + \frac{1}{2} g h^2 \end{bmatrix}_{,x} = \begin{bmatrix}[1.5]
\left( h_{\text{MS}} \right)_{,t} + \left( h_{\text{MS}} u_{\text{MS}} \right)_{,x} \\
\left(h_{\text{MS}} u_{\text{MS}} \right)_{,t} + \left( h_{\text{MS}} u^2_{\text{MS}} + \frac{1}{2} g h_{\text{MS}}^2\right)_{,x}\end{bmatrix}.
\end{equation}

We compute numerical solutions to \cref{eqn:forced_shllw} with \cref{eqn:meth_man_solns},
on $[-1,1]$, with periodic boundary conditions, with $\melems = N = 10 \times 2^k$ elements for $k=0,1,2,\ldots,5$,
to time $t=0.5$, and with three different orders of accuracy: $\morder=3,4,5$.
For each simulation, we aim to compute the relative  $L^2[-1,1]$ error:
\begin{equation}
\text{rel.} \, \,  L^2[-1,1] \, \text{error} := \sum_{\ell=1}^{\meq}
\sqrt{\frac{\int_{-1}^{1} \left| q^h_{\ell}(0.5,x) - q^{\star}_{\ell}(0.5,x) \right|^2 \, dx}{\int_{-1}^{1} \left| q^{\star}_{\ell}(0.5,x) \right|^2 \, dx}},
\end{equation}
where $q^h$ and $q^{\star}$ are the numerical and exact solutions, respectively.
In practice, we replace the exact solution by a piecewise Legendre polynomial approximation of
degree $\morder+1$, which allows us to obtain the following approximate relative $L^2[-1,1]$ error:
\begin{equation}
\label{eqn:rel_error}
e_N := \sum_{\ell=1}^{\meq} \sqrt{ \frac{\sum\limits_{i=1}^{N} \left( \sum\limits_{k=1}^{\mcorr} \left( Q_{i \, (k,\ell)} -
Q^{\star}_{i \, (k,\ell)} \right)^2 + \left( Q^{\star}_{i \, (\mcorr+1,\ell)} \right)^2 \right)}{
\sum\limits_{i=1}^{N} \sum\limits_{k=1}^{\mcorr+1} \left( Q^{\star}_{i \, (k,\ell)} \right)^2 } },
\end{equation}
where $Q^{\star}$ are the Legendre coefficients of the exact solution \cref{eqn:meth_man_solns}, written
as conservative variables at
time $t=0.5$. The exact solution coefficients are computed via Gaussian quadrature
with $\morder+1$ points:
\begin{equation}
\vec{Q^{\star}_{i \, (k,:)}} := \frac{1}{2} \sum_{a=1}^{\morder+1} \omega^{\star}_a \, \vec{q_{\text{MS}}}\left(0.5,x_i
+ \frac{\Delta x}{2} \mu^{\star}_a \right),
\end{equation}
where $\omega^{\star}_a$ and $\mu_a^{\star}$ for $a=1,\ldots,\morder+1$ are the weights and abscissas
of the $(\morder+1)$-point Gaussian quadrature rule. The resulting errors and error ratios are catalogued
in \cref{table:shllw_1d_error}.

\begin{table}
\begin{center}
\begin{Large}
\begin{tabular}{|c||c|c||c|c||c|c|}
\hline
{\normalsize $N$} & {\normalsize $\text{e}_{N} (\morder=3)$} & {\normalsize $\log_2\frac{\text{e}_{N/2}}{\text{e}_{N}}$}  & {\normalsize $\text{e}_N (\morder=4)$} & {\normalsize $\log_2\frac{\text{e}_{N/2}}{\text{e}_{N}}$}  & {\normalsize $\text{e}_N (\morder=5)$} & {\normalsize $\log_2\frac{\text{e}_{N/2}}{\text{e}_{N}}$} \\
\hline\hline
{\normalsize 10} & {\normalsize 1.990e-02} & -- & {\normalsize 1.324e-03} & -- & {\normalsize 1.104e-04} & -- \\\hline
{\normalsize 20} & {\normalsize 2.409e-03} & {\normalsize $3.047$} & {\normalsize 7.650e-05} & {\normalsize $4.114$} & {\normalsize 3.216e-06} & {\normalsize $5.102$} \\\hline
{\normalsize 40} & {\normalsize 3.183e-04} & {\normalsize $2.920$} & {\normalsize 4.497e-06} & {\normalsize $4.088$} & {\normalsize 1.090e-07} & {\normalsize $4.883$} \\\hline
{\normalsize 80} & {\normalsize 4.210e-05} & {\normalsize $2.918$} & {\normalsize 2.792e-07} & {\normalsize $4.010$} & {\normalsize 3.679e-09} & {\normalsize $4.889$} \\\hline
{\normalsize 160} & {\normalsize 5.563e-06} & {\normalsize $2.920$} & {\normalsize 1.768e-08} & {\normalsize $3.981$} & {\normalsize 1.199e-10} & {\normalsize $4.939$} \\\hline
{\normalsize 320} & {\normalsize 7.341e-07} & {\normalsize $2.922$} & {\normalsize 1.126e-09} & {\normalsize $3.972$} & {\normalsize 3.846e-12} & {\normalsize $4.963$} \\\hline
\end{tabular}
\caption{Relative $L^2$ errors for the forced shallow water equations with periodic boundary conditions.
\label{table:shllw_1d_error}}
\end{Large}
\end{center}
\end{table}

\subsubsection{Shock formation in finite time}
In \cref{subsec:burgers_example} we considered an example for the Burgers equation that started as a smooth initial condition and then shocked in finite time. We now attempt similar example, this time for the shallow
water equations. We take the following initial conditions on the computational domain $[-1,1]$ with outflow
boundary conditions at $x=\pm 1$:
\begin{equation}
\label{eqn:shllw_gauss_example}
h(t=0,x) = 1 + e^{-100x^2}, \quad u(t=0,x) = 0.
\end{equation}
We note that the subsequent dynamics will result in the initial Gaussian bump splitting into
two smaller height disturbances, one propagating to the left, the other to the right. The tops of these height disturbances will propagate faster than the rest of the profile, which will lead to wave steepening and eventually to
 the formation of a shock in both the left and right propagating disturbances.

The results of solving the shallow water equations with initial conditions \cref{eqn:shllw_gauss_example}, with the
$\morder=4$ version of the scheme, and with $\melems=100$ elements, is shown in \cref{fig:shllw_gaussian}.
 As in \cref{subsec:burgers_example}, we run this problem
with and without limiters in order to demonstrate efficacy of the limiting procedure.
The individual panels show at time $t=0.5$: (a) the height and (b) velocity as computed with the unlimited scheme,
and (c) the height and (d) velocity as computed with the limited scheme.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution. In both cases, the numerical solution is far from
exhibiting positivity violations; and therefore, the only limiter that is active is the
correction step limiter (see \cref{subsec:nonoscillatory}). We again see the ability of this once-per-time-step
limiter to simultaneously remove unphysical oscillations and to maintain the sharpness of the solution.

\begin{figure}[!t]
\centering
\begin{tabular}{cc}
(a)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_shllw_gaussian_height_no_limiter.pdf} &
(b)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_shllw_gaussian_velocity_no_limiter.pdf} \\
(c)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_shllw_gaussian_height.pdf} &
(d)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_shllw_gaussian_velocity.pdf}
\end{tabular}
\caption{Comparison of the unlimited and limited schemes on the shallow water equations \cref{eqn:shllw}.
Shown are solutions computed with an initially smooth height profile and zero velocity,
\cref{eqn:shllw_gauss_example}, that forms shocks in finite time, outflow boundary conditions,
$\melems = 100$, and $\morder=4$ (i.e., fourth-order scheme with polynomial degree $\mdeg=3$). The solutions are shown at time $t=0.5$, after the shocks have formed.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution. The individual panels show the (a) height from
the unlimited scheme, (b) velocity from the unlimited scheme, (c) height from the scheme with limiters,
and (d) velocity from the scheme with limiters.
 \label{fig:shllw_gaussian}}
\end{figure}

\subsubsection{Dambreak problem}
Next we consider an example of a Riemann problem, which for shallow water equations is also called the
dambreak problem. The initial conditions are piecewise constant:
\begin{equation}
\label{eqn:shllw_dambreak}
\left( h, \, u \right)(t=0,x) =
\begin{cases}
   \left( 1.0, \, 0 \right) & \quad x<0, \\
   \left( 0.1, \, 0 \right) & \quad x>0.
\end{cases}
\end{equation}
The entropy-satisfying solution with this initial data is a left-propagating rarefaction wave and a right-propagating
shock. A full mathematical explanation of the solution of this Riemann problem can be found in several textbooks, including Chapter 13 of LeVeque \cite{book:Le02}.

The results of solving the shallow water equations with initial conditions \cref{eqn:shllw_dambreak}, on
the domain $[-1,1]$ with outflow boundary conditions, with the
$\morder=4$ version of the scheme, and with $\melems=200$ elements is shown in \cref{fig:shllw_dambreak}.
The individual panels show at time $t=0.6$: (a) the height and (b) the velocity,
 with the exact Riemann solution superimposed.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution.
As in the previous example, the numerical solution is far from
exhibiting positivity violations; and therefore, the only limiter that is active is the
correction step limiter (see \cref{subsec:nonoscillatory}). We again see the ability of this once-per-time-step
limiter to simultaneously remove unphysical oscillations and to maintain the sharpness of the solution.

\begin{figure}[!t]
\centering
\begin{tabular}{cc}
(a)\includegraphics[width=0.43\textwidth,trim={0.7cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_shllw_dambreak_height.pdf} &
(b)\includegraphics[width=0.43\textwidth,trim={0.7cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_shllw_dambreak_velocity.pdf}
\end{tabular}
\caption{
Dambreak Riemann problem for the shallow water equations \cref{eqn:shllw}.
Shown are solutions computed with initial conditions
\cref{eqn:shllw_dambreak}, outflow boundary conditions,
$\melems = 200$, and $\morder=4$ (i.e., fourth-order scheme with polynomial degree $\mdeg=3$). The solutions are shown at time $t=0.4$.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution. The individual panels show the numerical solution and the superimposed exact Riemann solution for the (a) height and
(b) fluid velocity.
\label{fig:shllw_dambreak}}
\end{figure}

\subsubsection{Double rarefaction}
\label{subsubsec:shllw_double_raref}
In order to test the positivity limiters on the shallow water equations, we
attempt a different Riemann problem, this time with initial conditions
that result in two counter-propagating rarefactions that leave in their wake a near-vacuum
state. The initial conditions are as follows:
\begin{equation}
\label{eqn:shllw_double_raref}
\left( h, \, u  \right)(t=0,x) =
\begin{cases}
   \left( 1, \, -2 \right) & \quad x<0, \\
   \left( 1, \, +2 \right) & \quad x>0.
\end{cases}
\end{equation}
A full mathematical explanation of the solution of this Riemann problem can be found in several textbooks, including in Chapter 13 of LeVeque \cite{book:Le02}.

The results of solving the  shallow water equations with initial conditions \cref{eqn:shllw_double_raref}, on
the domain $[-1,1]$ with outflow boundary conditions, with the
$\morder=4$ version of the scheme, and with $\melems=200$ elements is shown in \cref{fig:euler_doublerarefaction}.
The individual panels show at time $t=0.25$: (a) the full height profile, (b) a zoomed-in view of the height in the
near-vacuum region, (c) the full momentum profile, and (d) a zoomed-in view of the
momentum. The exact Riemann solution is superimposed.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution.
The proposed limiters are able to handle the near-vacuum solution in the center of the computational domain,
and the numerical solution remains stable and positivity-preserving in the sense of \cref{eqn:shllw_pos_set_1}
and \cref{eqn:shllw_pos_set_2}.

\begin{figure}[!t]
\centering
\begin{tabular}{cc}
(a)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_shllw_doublerarefaction_height.pdf} &
(b)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_shllw_doublerarefaction_height_zoom.pdf} \\
(c)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_shllw_doublerarefaction_mom.pdf} &
(d)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_shllw_doublerarefaction_mom_zoom.pdf}
\end{tabular}
\caption{Double rarefaction Riemann problem for the shallow water equation \cref{eqn:shllw}.
Shown are solutions computed with initial conditions
\cref{eqn:shllw_double_raref}, outflow boundary conditions,
$\melems = 200$, and $\morder=4$ (i.e., fourth-order scheme with polynomial degree $\mdeg=3$).
The individual panels show at time $t=0.25$: (a) the full height profile, (b) a zoomed-in view of the height in the
near-vacuum region, (c) the full momentum profile, and (d) a zoomed-in view of the
momentum. The exact Riemann solution is superimposed.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution.
\label{fig:shllw_doublerarefaction}}
\end{figure}

\subsection{\it Compressible Euler equations}
Finally we consider solving the compressible Euler equations \cref{eqn:euler} with various initial conditions that
demonstrate the accuracy of the scheme and the efficacy of the proposed limiting strategies.
In all examples the adiabatic constant is taken as $\gamma=1.4$ and the limiter parameter is taken as
$\varepsilon = 10^{-14}$.

\subsubsection{Convergence test}
The compressible Euler equations, along with suitable boundary conditions,
admit a class of non-trivial and non-stationary
exact solutions in which the fluid velocity and pressure remain globally constant, but
an arbitrary density profile is advected by the fluid velocity.
We will use such a solution to verify the order of accuracy of the proposed scheme.
In particular, we consider the following exact solution on $[-1,1]$ with periodic boundary conditions:
\begin{equation}
\label{eqn:euler_smooth_ic}
\rho(t,x) = 1 + 0.5 \sin\left(3 \pi \left( x - 0.5 t \right) \right), \quad
u(t,x) = 0.5, \quad p(t,x) = 0.75.
\end{equation}

We compute numerical solutions to \cref{eqn:euler} with initial conditions
obtained from \cref{eqn:euler_smooth_ic},
on $[-1,1]$, with periodic boundary conditions, with $\melems = N = 10 \times 2^k$ elements for $k=0,1,2,\ldots,5$, to time $t=1$, and with three different orders of accuracy: $\morder=3,4,5$. For each simulation we compute
the relative $L^2[-1,1]$ error via \cref{eqn:rel_error}, where the $Q^{\star}$ coefficients
are computed from the conservative variable version of \cref{eqn:euler_smooth_ic} at time $t=1$.
 The resulting errors and error ratios are catalogued
in \cref{table:euler_1d_error}.

\begin{table}
\begin{center}
\begin{Large}
\begin{tabular}{|c||c|c||c|c||c|c|}
\hline
{\normalsize $N$} & {\normalsize $\text{e}_N (\morder=3)$} & {\normalsize $\log_2\frac{\text{e}_{N/2}}{\text{e}_{N}}$}  & {\normalsize $\text{e}_N (\morder=4)$} & {\normalsize $\log_2\frac{\text{e}_{N/2}}{\text{e}_{N}}$}  & {\normalsize $\text{e}_N (\morder=5)$} & {\normalsize $\log_2\frac{\text{e}_{N/2}}{\text{e}_{N}}$} \\
\hline\hline
{\normalsize 10} & {\normalsize 2.161e-02} & -- & {\normalsize 3.109e-03} & -- & {\normalsize 2.179e-04} & -- \\\hline
{\normalsize 20} & {\normalsize 3.742e-03} & {\normalsize $2.530$} & {\normalsize 1.225e-04} & {\normalsize $4.665$} & {\normalsize 1.010e-05} & {\normalsize $4.431$} \\\hline
{\normalsize 40} & {\normalsize 6.540e-04} & {\normalsize $2.517$} & {\normalsize 7.182e-06} & {\normalsize $4.093$} & {\normalsize 4.438e-07} & {\normalsize $4.509$} \\\hline
{\normalsize 80} & {\normalsize 9.633e-05} & {\normalsize $2.763$} & {\normalsize 4.398e-07} & {\normalsize $4.029$} & {\normalsize 1.623e-08} & {\normalsize $4.773$} \\\hline
{\normalsize 160} & {\normalsize 1.279e-05} & {\normalsize $2.913$} & {\normalsize 2.728e-08} & {\normalsize $4.011$} & {\normalsize 5.343e-10} & {\normalsize $4.925$} \\\hline
{\normalsize 320} & {\normalsize 1.629e-06} & {\normalsize $2.973$} & {\normalsize 1.706e-09} & {\normalsize $3.999$} & {\normalsize 1.695e-11} & {\normalsize $4.979$} \\\hline
\end{tabular}
\caption{Relative $L^2$ errors for the compressible Euler equations with constant pressure and fluid velocity and with periodic boundary conditions.
\label{table:euler_1d_error}}
\end{Large}
\end{center}
\end{table}

\subsubsection{Shock tube problem}
Next we consider an example of a Riemann problem, which for compressible Euler equations is also called the
shock tube problem. We consider the celebrated Sod shock tube problem \cite{article:Sod78}, for
which the initial conditions are
\begin{equation}
\label{eqn:euler_shocktube}
\left( \rho, \, u, \, p \right)(t=0,x) =
\begin{cases}
   \left( 1.000, \, 0, \, 1.0 \right) & \quad x<0, \\
   \left( 0.125, \, 0, \, 0.1 \right) & \quad x>0.
\end{cases}
\end{equation}
The entropy-satisfying solution with this initial data is a left-propagating rarefaction wave,
a right-propagating contact wave, and a (faster) right-propagating
shock. A full mathematical explanation of the solution of this Riemann problem can be found in several textbooks, including Chapter 14 of LeVeque \cite{book:Le02}.

The results of solving the shallow water equations with initial conditions \cref{eqn:euler_shocktube}, on
the domain $[-1,1]$ with outflow boundary conditions, with the
$\morder=4$ version of the scheme, and with $\melems=200$ elements is shown in \cref{fig:euler_shocktube}.
The individual panels show at time $t=0.4$: (a) the density, (b) velocity, and (b) the pressure,
 with the exact Riemann solution superimposed.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution.
The numerical solution is far from
exhibiting positivity violations; and therefore, the only limiter that is active is the
correction step limiter (see \cref{subsec:nonoscillatory}). This example
demonstrates the ability of this once-per-time-step
limiter to simultaneously avoid producing unphysical oscillations and
to maintain the sharpness of the solution.

\begin{figure}[!t]
\centering
\begin{tabular}{cc}
(a)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_euler_shocktube_density.pdf} &
(b)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_euler_shocktube_velocity.pdf} \\
\multicolumn{2}{c}{(c)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_euler_shocktube_pressure.pdf}}
\end{tabular}
\caption{
Sod shocktube Riemann problem for the compressible Euler equations \cref{eqn:euler}.
Shown are solutions computed with initial conditions
\cref{eqn:euler_shocktube}, outflow boundary conditions,
$\melems = 200$, and $\morder=4$ (i.e., fourth-order scheme with polynomial degree $\mdeg=3$). The solutions are shown at time $t=0.4$.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution. The individual panels show the numerical solution and the superimposed exact Riemann solution for the (a) density,
(b) fluid velocity, and (c) pressure.
\label{fig:euler_shocktube}}
\end{figure}

\subsubsection{Double rarefaction}
In order to test the positivity limiters on the compressible Euler equations, we again
attempt a Riemann problem that results in two counter-propagating rarefactions that leave in their wake a near-vacuum
state (see \cref{subsubsec:shllw_double_raref} for the shallow water version).
This particular test case can be found in several papers, e.g., \cite{article:ZhangShu10rectmesh,article:ZhangShu11,article:seal2014explicit,article:MoRoSe17}.
The initial conditions are as follows:
\begin{equation}
\label{eqn:euler_double_raref}
\left( \rho, \, u, \, p \right)(t=0,x) =
\begin{cases}
   \left( 7, \, -1, \, 0.2 \right) & \quad x<0, \\
   \left( 7, \, +1, \, 0.2 \right) & \quad x>0.
\end{cases}
\end{equation}
A full mathematical explanation of the solution of this Riemann problem can be found in several textbooks, including Chapter 14 of LeVeque.\cite{book:Le02}.

The results of solving the compressible Euler equations with initial conditions \cref{eqn:euler_double_raref}, on
the domain $[-1,1]$ with outflow boundary conditions, with the
$\morder=4$ version of the scheme, and with $\melems=200$ elements is shown in \cref{fig:euler_doublerarefaction}.
The individual panels show the (a) density, (b) velocity, and (c) pressure at time $t=0.6$ with the exact
Riemann solution superimposed.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution.
The proposed limiters are able to handle the near-vacuum solution in the center of the computational domain,
and the numerical solution remains stable and positivity-preserving in the sense of \cref{eqn:euler_pos_set_1}
and \cref{eqn:euler_pos_set_2}.

\begin{figure}[!t]
\centering
\begin{tabular}{cc}
(a)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_euler_doublerarefaction_density.pdf} &
(b)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_euler_doublerarefaction_velocity.pdf} \\
\multicolumn{2}{c}{
(c)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_euler_doublerarefaction_pressure.pdf}}
\end{tabular}
\caption{
Double rarefaction Riemann problem for the compressible Euler equations \cref{eqn:euler}.
Shown are solutions computed with initial conditions
\cref{eqn:euler_double_raref}, outflow boundary conditions,
$\melems = 200$, and $\morder=4$ (i.e., fourth-order scheme with polynomial degree $\mdeg=3$). The solutions are shown at time $t=0.6$.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution. The individual panels show the numerical solution and the superimposed exact Riemann solution for the (a) density,
(b) fluid velocity, and (c) pressure.
\label{fig:euler_doublerarefaction}}
\end{figure}

\subsubsection{Sedov blast problem}
Another standard test case for verifying the efficacy of positivity-preserving
limiters for the Euler equations is the 1D Sedov blast problem \cite{book:sedov1993similarity}.
For example, this test case is featured in the following papers:
\cite{article:ZhangShu10rectmesh,article:MoRoSe17,article:ZhangShu12}.
The initial conditions for the problem can be written on $[-1,1]$ as
\begin{equation}
\label{eqn:sedov_ic}
\rho(t=0,x) = 1, \quad
u(t=0,x) = 0, \quad
p(t=0,x) = \begin{cases}
(\gamma-1)\left(\frac{3.2 \times 10^6}{\Delta x} \right) & \bigl| x \bigr| \le \frac{\Delta x}{2}, \\
(\gamma-1) \left(10^{-12} \right) & \text{otherwise},
\end{cases}
\end{equation}
where $\Delta x$ is the mesh grid spacing. We assume here that there are an odd number of mesh elements
so that the middle element, $i=(\melems+1)/2$, is centered at the origin and defined by
\begin{equation}
{\mathcal T}_{(\melems+1)/2} = \left[-\frac{\Delta x}{2}, \frac{\Delta x}{2} \right].
\end{equation}
The initial conditions \cref{eqn:sedov_ic} represent a constant solution with almost zero
pressure everywhere, except in the middle element, where the pressure is many orders of magnitude larger than its surroundings. The resulting solution is a pressure blast that emenates
from the central element and propagates outward in both directions.
These
initial conditions approximate a delta function of pressure.  Once the wave front propagates
away from the center of the domain, a post-shock region with near zero density is left behind.
This example can only be simulated with methods that either use extremely aggressive limiters or methods with guaranteed positivity-preservation.

The results of solving the compressible Euler equations with initial conditions \cref{eqn:sedov_ic}, on
the domain $[-1,1]$ with outflow boundary conditions, with the
$\morder=4$ version of the scheme, and with $\melems=201$ elements is shown in \cref{fig:euler_sedov}.
The individual panels show the (a) density, (b) velocity, and (c) pressure at time $t=0.0004$.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution.
The proposed limiters are able to handle the massive pressure jumps and the
near-vacuum solution in the center of the computational domain.
The numerical solution remains stable and positivity-preserving in the sense of \cref{eqn:euler_pos_set_1}
and \cref{eqn:euler_pos_set_2}.

\begin{figure}[!t]
\centering
\begin{tabular}{cc}
(a)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_euler_sedov_density.pdf} &
(b)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_euler_sedov_velocity.pdf} \\
\multicolumn{2}{c}{(c)\includegraphics[width=0.43\textwidth,trim={1.0cm 1.8cm 1.3cm 1.8cm},clip]{figures/num_euler_sedov_pressure.pdf}}
\end{tabular}
\caption{1D Sedov blast wave problem for the compressible Euler equations \cref{eqn:euler}.
Shown are solutions computed with initial conditions
\cref{eqn:sedov_ic}, outflow boundary conditions,
$\melems = 201$, and $\morder=4$ (i.e., fourth-order scheme with polynomial degree $\mdeg=3$). The solutions are shown at time $t=0.0004$.
In each panel we are plotting four points per element in order to clearly show
the subcell structure of the numerical solution. The individual panels show the (a) density,
(b) fluid velocity, and (c) pressure.
\label{fig:euler_sedov}}
\end{figure}


TITLE: $\cos^2 \alpha + \cos^2 \beta + \cos^2 \gamma =1$
QUESTION [2 upvotes]: Let be $\alpha, \beta, \gamma$  the angles between a generic direction in 3D and the axes $x,y,z$, respectively.
Prove that
$\cos^2 \alpha + \cos^2 \beta + \cos^2 \gamma =1$.

PS: the 2D case is trivial. But I can't prove the 3D case.

REPLY [3 votes]: Consider the unit vector $$(u,v,w).$$
The cosines of the angles it forms with the axis are given by the dot products
$$\cos\alpha=(1,0,0)\cdot(u,v,w)=u,\\\cos\beta=(0,1,0)\cdot(u,v,w)=v,\\\cos\gamma=(0,0,1)\cdot(u,v,w)=w,$$
and as it is unit
$$u^2+v^2+w^2=1.$$


TITLE: Rhombus tilings with more than three directions
QUESTION [20 upvotes]: The point of this question is to construct a list of references on the following subject: Fix vectors $v_1$, $v_2$, ..., $v_g$ in $\mathbb{R}^2$, all lying in a half plane in that cyclic order. I am interested in tilings by parallelograms of the form $\mathrm{Hull}(z, z+v_i, z+v_i+v_j, z+v_j)$. It is traditional in this subject to take all the $v_i$ of the same length, so that the tiles are rhombi. In my personal opinion, this convention is unmotivated, but I will obey convention and speak of rhombus tilings in the rest of this post.
The most obvious region to tile is a centrally symmetric $2g$-gon, with edges of the form $a_1 v_1$, $a_2 v_2$, ..., $a_g v_g$, $-a_1 v_1$, ..., $- a_g v_g$ for some positive integers $a_i$. One might also study tilings of the whole plane, perhaps with some periodicity conditions. I'm also open to hearing about results on regions of other shapes.
What I do want to do is rule out results which are special to the case $g=3$. Rhombus tilings
of a centrally symmetric hexagon are in bijection with plane partitions, and hence all the tools of the plane partition literature are available. They are also related to perfect matchings of the hexagonal grid, so they can be studied by Kastelyn's method and related tools. They are also in bijection with certain families of non-crossing paths, so they can be studied by the Gessel-Viennot method. There is an embarrassment of tools available for the $g=3$ case. I want to limit this discussion only to tools that work for higher $g$.
Below the line is what I know so far:

$\bullet$ The two most obvious special cases are (1) tilings of an octagon where $(a_1, a_2, a_3, a_4) = (k,k,k,k)$ and (2) tilings of a $2g$-gon with $a_1=a_2=\cdots=a_g=1$. These are counted by Sloane A093937 and A006245 respectively. There is no closed formula in either case.
$\bullet$ Such rhombus tilings are in bijection with reduced words in $S_{a_1+a_2+\cdots+a_g}$ modulo commutation. In particular, the case where $a_1=a_2=\cdots=a_g=1$ corresponds to reduced words for the long word in $S_g$, modulo commutation. This has been rediscovered several times, but the canonical reference seems to be Rhombic tilings of polygons and classes of reduced words in Coxeter groups, by Serge Elnitsky.
$\bullet$ Drawing the paths dual to the rhombi gives us an arrangement of pseudo-lines, sometimes called deBruijn paths.
$\bullet$ Pseudo-line arrangements can be described by oriented matroids. Tracing through this language, the $a_1=a_2=\cdots=a_g=1$ case corresponds to rank $3$ oriented matroids $M$ on the set $\{ 1,2,\ldots, g, \infty \}$ such that the underlying matroid is uniform and $M/\infty$ is the uniform rank $2$ oriented matroid where $\{1,2,\ldots,g \}$ appear in that cyclic order. Call that rank $2$ matroid $U(g)$. Dualizing, we are looking for extensions of the rank $g-2$ oriented matroid $U(g)^{\ast}$ whose underlying matroid is uniform. One can generalize this description to the general case: The matroid $U(g)$ becomes the oriented matroid of rank $2$ made up of $g$ parallelism classes, of sizes $a_1$, ..., $a_g$, in that cyclic order.
Extensions of oriented matroids were studied by Sturmfels and Ziegler. In particular, they proved that a certain poset associated to this extension problem is homotopic to a sphere.
$\bullet$ Let $A$ and $B$ be two subsets of $\{ 1,2,\ldots, g \}$. We say that $A$ and $B$ are strongly separated if either every element of $A \setminus B$ is less than every element of $B \setminus A$, or is every element of $B \setminus A$ is less than every element of $A \setminus B$. A collection $\mathcal{C}$ of subsets of $\{ 1,2,\ldots, g \}$ is called strongly separated if every two elements in the collection are strongly separated. A collection is called maximally strongly separated if it is strongly separated and is not contained in any larger strongly separated collection.
Maximal strongly separated sets are in bijection with rhombus tilings of the $(1,1,\ldots,1)$-gon by a result of Leclerc and Zelevinsky. Although it doesn't appear to be in the literature, the generalization to other centrally symmetric polygons is easy enough: One uses the multiset with $a_i$ copies of $i$. Here if $A$ and $B$ are multisets which contain $u$ and $v$ copies of $i$, then $A \setminus B$ is the multiset which contains $\max(u-v, 0)$ copies of $i$.
$\bullet$ Any two rhombus tilings of the same $2g$-gon are linked by a sequence of flips, where a flip takes a single hexagon and retiles it. This is the same as the fact that any two reduced words are (up to commutation) joined by a sequence of braid moves.
$\bullet$ Take the flip graph referred to in the previous bullet point and direct it so that, in the language of braid moves, an edge points from $\cdots s_i s_{i+1} s_i \cdots$ to $\cdots s_{i+1} s_i s_{i+1} \cdots$. The resulting graph is acyclic, and is the Hasse diagram of the poset obtained by taking the transitive closure.
This poset is graded and has unique minimum and maximum elements.  When $a_1=a_2=\cdots=a_g=1$, it is called the higher Bruhat order $B(g,2)$. It is a lattice for small $g$, but $B(6,2)$ is not a lattice.
$\bullet$ Let $1 \leq i < j < k \leq n$ and let $\gamma_i$, $\gamma_j$ and $\gamma_k$ be pseudolines in directions $(i,j,k)$. Then they either intersect each other in the order $(i,j)$, $(i,k)$, $(j,k)$, or the reverse of this order. A triple of such pseudo-lines is called a triangle. The triangle is ordinary if the intersections occur in the first manner, and inverted if the occur in the second.
For $\gamma$ and $\delta$ in $B(g,2)$, we have $\gamma \leq \delta$ if and only if every triangle which is inverted in $\gamma$ is also inverted in $\delta$. Presumably, this is also true when some of the sides of the polygon are longer than $1$, but I don't know a reference.
For $\gamma$ and $\delta$ two pseudo-line arrangements, the flip distance from $\gamma$ to $\delta$ is bounded below by the number of triangles which are inverted in one of $\gamma$ and $\delta$ but not the other. This lower bound is exact for $g=3$ or $4$, for $g=5$ if one of the $a_i$ is $1$, and for $(2,2,2,2,2)$, but not for any other case.
$\bullet$ The Bohne-Dress theorem: Let $C$ be the cube $\{ (x_1, \ldots, x_g) : 0 \leq x_i \leq a_i \} \subset \mathbb{R}^g$. Map $C$ to $\mathbb{R}^2$ by the linear map $\pi$ sending the $i$-th basis vector to $v_i$. So the $2g$-gon is the image of $\pi$. Then rhombus tilings of $\pi(C)$ are in bijection with continuous sections $\sigma: \pi(C) \to C$ whose image lies in the $2$-skeleton of the obvious cubical complex.
$\bullet$ Miscellaneous: some interesting lemmas about rhombus tilings are proved in my paper with Andre Henriques and in Shapiro-Shapiro-Vainshtein.

Some things I particularly don't know: Do we know the approximate number of these tilings? Can we sample uniformly from this set? Has anyone worked on periodic rhombus tilings? (Surely someone has, but I haven't found them.) What do random tilings of a large octagon look like? (Even if nothing is proven, someone presumably has a guess.) Is there a way of comparing elements of $B(g,2)$ in time less than order of $g^3$?

REPLY [8 votes]: This post is almost ten years old, but I find it interesting and I would like to add a few pointers from my PhD, itself about 20 years old.
First, in the paper Integer Partitions, Tilings of 2D-gons and Lattices
 we explored the relations between rhombus tilings and integer partitions, and obtained a decomposition into disjoint distributive lattices of the flip graph of any 2D-gon tilings. In addition, the quotient of this set by these lattices has itself the same structure. These are direct consequences of the fact that they are isomorphic to generalized integer partitions.
Here is a simple example for the flip graph of an elemenraty octogon:

Second, in the paper Graph encoding of 2D-gon tilings, we used the tile adjacency graph to process rhombus tilings and sample them, like the random decagon tiling below (the full resolution vector file is available on demand):


TITLE: Differentiate $y =\sin(1+x^2)^{1/2}$
QUESTION [0 upvotes]: I've tried differentiating $y= \sin(1+x^2)^{1/2}$ using the chain rule, but I keep getting the wrong answer. Can anyone give me a step by step so I can see what I'm doing wrong? Thanks.

REPLY [0 votes]: If $y=\sqrt{u}$, $u=\sin(v)$, and $v=1+x^2$, then
$$\frac{dy}{dx}=\frac{dy}{du}\cdot \frac{du}{dv}\cdot \frac{dv}{dx}$$
So,
$$\frac{dy}{dx}=\frac{1}{2\sqrt{u}}\cdot \cos(v)\cdot(2x)$$
Expressing this in terms of $x$, we get:
$$\frac{dy}{dx}=\frac{1}{2\sqrt{\sin(1+x^2)}}\cdot \cos(1+x^2)\cdot(2x)=\frac{x\cos(1+x^2)}{\sqrt{\sin(1+x^2)}}$$
This is similar with the answers of the guys previously, where the idea is the "Chain Rule". I have posted this answer because one may find it more clear.