text:
\begin{document}

\title{A high-order relaxation method with projective integration for solving nonlinear systems of hyperbolic conservation laws}

\author{
Pauline Lafitte
	\thanks{Laboratoire de Math\'ematiques Appliqu\'ees aux Syst\`emes, Ecole Centrale Paris, Grande Voie des Vignes, 92290 Ch\^atenay-Malabry, France ({\tt pauline.lafitte@ecp.fr}).} \and
Ward Melis
	\thanks{Department of Computer Science, K.U. Leuven, Celestijnenlaan 200A, 3001 Leuven, Belgium ({\tt ward.melis@cs.kuleuven.be}).} \and
Giovanni Samaey
	\thanks{Department of Computer Science, K.U. Leuven, Celestijnenlaan 200A, 3001 Leuven, Belgium ({\tt giovanni.samaey@cs.kuleuven.be}).}}

\maketitle

\begin{abstract}
We present a general, high-order, fully explicit relaxation scheme for systems of nonlinear hyperbolic conservation laws in multiple dimensions. The scheme consists of two steps: first, the nonlinear hyperbolic conservation law is approximated by a kinetic equation with stiff BGK source term. Then, this kinetic equation is integrated in time using a projective integration method. After taking a few small (inner) steps with a simple, explicit method (such as direct forward Euler) to damp out the stiff components of the solution, the time derivative is estimated and used in an (outer) Runge-Kutta method of arbitrary order.  We show that, with an appropriate choice of inner step size, the time step restriction on the outer time step is similar to the CFL condition for the hyperbolic conservation law. Moreover, the number of inner time steps is also independent of the stiffness of the BGK source term. We discuss stability and consistency, and illustrate with numerical results (linear advection, Burgers' equation and the shallow water and Euler equations) in one and two spatial dimensions.
\end{abstract}

\section{Introduction} \label{sec:intro}
Hyperbolic conservation laws arise in numerous physical applications, such as fluid dynamics, plasma physics, traffic modeling and electromagnetism (see, e.g., \cite{whitham2011linear,leveque2002finite}). They express the conservation of physical quantities (such as mass, momentum, or energy) and may be supplemented with boundary conditions that control influx or outflux at the boundaries of the physical domain \cite{leveque2002finite}. In this paper, we consider a system of hyperbolic conservation laws in multiple spatial dimensions:
\begin{equation} \label{eq:cons_law_diff_system} \partial_t \u + \nabla_{\x} \cdot \F(\u)= 0, \end{equation}
or, equivalently,
\begin{equation} \label{eq:cons_law_diff_system_long} \partial_t \u + \sum_{d=1}^{D}\partial_{x^d} \F^d(\u)= 0, \end{equation}
in which $\x=(x^d)_{d=1}^D \in \R^D$ represents the space variables ($D$ being the number of spatial dimensions), ${\u}(\x,t):=  (u_m(\x,t))_{m\in\{1,...,M\}} \in \R^{M}$ denotes the conserved quantities, and $\F(\u) \in \R^{M\times D}$ corresponds to the flux functions.

Hyperbolic conservation laws are often solved using a finite volume method \cite{leveque2002finite,morton2005numerical}, which is derived from the integral expression of the conservation law. To this end, in a scalar one-dimensional setting and with a spatially uniform grid, the domain is divided in $I$ cells $\mathcal{C}_i = [x_{i-1/2},x_{i+1/2}]$ with constant cell width $\Delta x$ over which the cell average of the solution $u(x,t)$ to the conservation law
\begin{equation} \label{eq:cons_law_diff_scalar} \partial_t u + \partial_x F(u)= 0, \end{equation}
is approximated at time $t=t^n$ by
\begin{equation} \label{fin_diff_U} U_i^n \approx \frac{1}{\Delta x}\int_{\mathcal{C}_i} u(x,t^n)dx. \end{equation}
Note that boldface is removed whenever the quantities are scalar. A numerical scheme is then constructed by integrating the conservation law \eqref{eq:cons_law_diff_scalar} in space over the cell $\mathcal{C}_i$ and in time from $t^n$ to $t^{n+1}$ to obtain
\begin{equation} \label{fin_diff_sch} U_i^{n+1} = U_i^n - \frac{\Delta t}{\Delta x}\left(F_{i+1/2}^n - F_{i-1/2}^n\right), \end{equation}
in which $\Delta t = t^{n+1} - t^n$ and the \emph{numerical flux} satisfies
\begin{equation} \label{fin_diff_F} F_{i\pm1/2}^n \approx \frac{1}{\Delta t}\int_{t^n}^{t^{n+1}} F\left(u\left(x_{i\pm1/2},t\right)\right)dt. \end{equation}
Clearly, equation \eqref{fin_diff_sch} is conservative by construction. The numerical fluxes $F_{i\pm1/2}^n$ can be obtained by constructing an (approximate) Riemann solver, based on a (possibly high-order) reconstruction of the solution in each of the cells using interpolation over the neighboring cells \cite{leveque2002finite,eno_weno}.

To avoid the (possibly tedious) computation of the solutions of local Riemann problems, relaxation methods offer an interesting alternative, see, e.g., discrete kinetic schemes in \cite{liu1987hyperbolic, jin1998diffusive, jin1995relaxation} and, in particular, \cite{aregba2000discrete} which also contains a brief historical overview. In a relaxation method, the conservation law \eqref{eq:cons_law_diff_system} is approximated by a problem of higher dimension containing a small relaxation parameter $\epsi$ such that, when $\epsi$ tends to zero, the original problem is recovered. The idea is that some of the difficulties associated with the original problem are avoided, while, for sufficiently small $\epsi$, the relaxation problem is a good approximation of the problem of interest. In this paper, we will consider the relaxation problem to be a kinetic BGK equation. In a scalar one-dimensional setting, this equation describes the evolution of a distribution function $f(x,v,t)$ of particles at position $x$ with velocity $v$ at time $t$ and takes the following form:
\begin{equation} \label{eq:kin_eq_1d_scalar} \partial_t f^\epsi + v\partial_x f^\epsi = \frac{1}{\epsi}\left(\Mv(u^\epsi) - f^\epsi\right). \end{equation}
The left hand side of equation \eqref{eq:kin_eq_1d_scalar} describes the transport of the particles whereas the right hand side represents the collisions between particles, which is modeled as a linear relaxation to the Maxwellian $\Mv(u^\epsi)$ with a relaxation time $\epsi$. The advantage of the kinetic equation \eqref{eq:kin_eq_1d_scalar} over the conservation law \eqref{eq:cons_law_diff_scalar} is the fact that the advection term in \eqref{eq:kin_eq_1d_scalar} is now linear; the disadvantage is the appearance of a stiff source term, which requires special care during time integration. The first methods, proposed in \cite{jin1995relaxation,aregba2000discrete} are based on splitting, thus restricting the order in time to 2. More recently, several \emph{asymptotic-preserving} methods based on IMEX techniques (in the sense of Jin \cite{jin1999efficient}) have been proposed that integrate the Boltzmann equation in the hyperbolic and diffusive regimes with a computational cost that is independent of $\epsi$ (see \cite{filbet2010class} and references within).  An appealing idea along this line of thought, based on IMEX Runge-Kutta methods, is presented in \cite{boscarino2013implicit}. Unfortunately, the proposed method is not very robust since it breaks down in intermediate regimes. An improvement was proposed in \cite{dimarco2013asymptotic}.

A robust and fully explicit alternative to splitting and IMEX, which allows for time integration of stiff systems with arbitrary order of accuracy in time is projective integration. Projective integration was proposed in \cite{gear2003projective} for stiff systems of ordinary differential equations. In such stiff problems, the fast modes, corresponding to the Jacobian eigenvalues with large negative real parts, decay quickly, whereas the slow modes correspond to eigenvalues of smaller magnitude and are the solution components of practical interest. Projective integration allows a stable yet explicit integration of such problems by first taking a few small (inner) steps with a simple, explicit method, until the transients corresponding to the fast modes have died out, and subsequently projecting (extrapolating) the solution forward in time over a large (outer) time step. In \cite{lafitte2012asymptotic}, projective integration was analyzed for kinetic equations with a diffusive scaling. An arbitrary order version, based on Runge-Kutta methods, has been proposed recently in \cite{lejon13}, where it was also analyzed for kinetic equations with an advection-diffusion limit. These methods fit within recent research efforts on numerical methods for multiscale simulation \cite{EEng03,E:2007p3747,KevrGearHymKevrRunTheo03,Kevrekidis2009a}; see also \cite{logg,sommeijer1990increasing,Vandekerckhove2007a} for related approaches. Alternative approaches to obtain a higher-order projective integration scheme have been proposed in \cite{Lee2007,Rico-Martinez}.

In this paper, we construct a relaxation method with projective integration to simulate arbitrary system of hyperbolic conservation laws in multiple space dimensions. The resulting scheme turns out to be fully explicit, of arbitrary order in time and fully general, avoiding the construction of complicated and problem specific (approximate) Riemann solvers. High-order projective Runge-Kutta methods are proposed and analyzed in \cite{lejon13}. The present paper focuses on their use in combination with a relaxation method, including ample numerical illustrations.

The remainder of this paper is structured as follows. In section \ref{sec:kinetic}, we will introduce the kinetic equations that will form the basis of the relaxation method, and discuss their asymptotic equivalence with the original hyperbolic problem. In section \ref{sec:proj_int}, we describe the projective integration method that will be used to integrate these kinetic equations. Based on the analysis in \cite{lejon13}, we then discuss convergence of the resulting projective integration relaxation method in the specific setting of hyperbolic conservation laws in section \ref{sec:stability}, including the choice of appropriate method parameters. Section \ref{sec:results} reports the results of extensive numerical tests for a set of benchmark problems in both one and two space dimensions: linear advection, nonlinear conservation, the dam-break problem and Sod's shock test. We conclude in section \ref{sec:conclusions} with a brief discussion and some ideas for future work.

\section{Relaxation systems} \label{sec:kinetic}

\subsection{Kinetic equation and hydrodynamic limit}
To solve equation \eqref{eq:cons_law_diff_system}, we introduce, as in \cite{aregba2000discrete}, the (hyperbolically scaled) kinetic equation
\begin{equation} \label{eq:kin_eq_short} \partial_t\f^\epsi + \v\cdot\nabla_{\x}\cdot \f^\epsi = \frac{1}{\epsi}\left(\Maxv(\u^\epsi) - \f^\epsi\right), \end{equation}
or, equivalently,
\begin{equation} \label{eq:kin_eq} \partial_t\f^\epsi + \sum_{d=1}^{D} v^d\partial_{x^d}\f^\epsi= \frac{1}{\epsi}\left(\Maxv(\u^\epsi) - \f^\epsi\right), \end{equation}
modeling the evolution of a vector of particle distribution functions $\f^\epsi(\x,\v,t)=(f^\epsi_m(\x,\v,t))_{m=1}^M \in~\R^M$. The particle positions and velocities are represented as $\x=(x^d)_{d=1}^D \in \R^D$ and $\v=(v^d)_{d=1}^D \in V \subset \R^D$ respectively, and the right hand side of \eqref{eq:kin_eq} represents a BGK collision operator \cite{bhatnagar1954model}, modeling linear relaxation of $\f^\epsi$ to a Maxwellian distribution $\Maxv(\u^\epsi)\in\R^M$, in which the argument $\u^\epsi(\x,t)=\langle\f^\epsi(\x,\v,t)\rangle$ is the density, obtained via averaging over the measured velocity space $(V,\mu)$,
\begin{equation}\label{eq:bracket} \u := \langle\f\rangle = \int_V \f d\mu(\v). \end{equation}

The advantage of this kinetic formulation is that the advection term is now linear, and therefore easier to discretize. The disadvantage is the increased dimension, as well as the introduction of the stiff source term of size $O(1/\epsi)$. The projective integration scheme that we will propose in section \ref{sec:proj_int} allows to integrate this stiff source term using an explicit method of arbitrary order.

To ensure that the kinetic equation \eqref{eq:kin_eq} converges to the conservation law \eqref{eq:cons_law_diff_system} in the hydrodynamic limit $\epsi \rightarrow 0$, one requires
\begin{equation}\label{eq:maxwellian_conditions}
  \begin{cases}
    \langle\Maxv(\u)\rangle = \u, \\
    \langle v^d \Maxv(\u)\rangle = \F^d(\u), \qquad 1 \le d \le D.
  \end{cases}
\end{equation}
Then, one can show \cite{aregba2000discrete} that, in the limit of $\epsi \to 0$, the kinetic model \eqref{eq:kin_eq} is approximated by the following equation:
\begin{equation} \label{eq:cons_law_limit_short} \partial_t\u^\epsi + \nabla_{\x}\cdot\F(\u^\epsi) = \epsi\; \nabla_{\x}\cdot\left(\B\;\nabla_{\x}\u^\epsi\right), \end{equation}
or, equivalently,
\begin{equation} \label{eq:cons_law_limit} \partial_t\u^\epsi + \sum_{d=1}^{D}\partial_{x^d} \F^d(\u^{\epsi})= \epsi\; \sum_{d=1}^D\partial_{x^d} \left(\sum_{d'=1}^D\B_{dd'}\;\partial_{x^{d'}}\u^\epsi\right), \end{equation}
with the diffusion matrix $\B$ given as
\begin{equation} \B_{dd'}(\u):=\langle v^dv^{d'}\partial_{\u}\Maxv(\u)\rangle - \partial_{\u}\F^{d} \partial_\u\F^{d'}, \end{equation}
in which the $M\times M$ matrices $\partial_{\u}\Maxv(\u)$ and $\partial_{\u}\F^{d}$ represent the Jacobian matrices of $\Maxv(\u)$ and $\F(\u)$ respectively.

Clearly, equation \eqref{eq:kin_eq} is consistent with equation \eqref{eq:cons_law_diff_system} to order $1$ in $\epsi$. Moreover, the analysis reveals an additional condition on $\Maxv$ and $V$. Indeed, to ensure the parabolicity of \eqref{eq:cons_law_limit}, the diffusion matrix $\B$ should be positive definite. This leads to the so-called \emph{subcharacteristic condition} \cite{bouchut1999construction,aregba2000discrete},
\begin{equation}\label{eq:subchar} \sum_{d,d'=1}^D\left(\B_{dd'}(\u)\xi^{d'}\cdot\xi^d\right)\ge 0, \end{equation}
for all $\xi^d$, $1\le d \le D$ in $\R^{M}$.

In what follows, we will always assume that the velocity space is discrete and of the form
\begin{equation} V := \left\{\v_{j}\right\}_{j=1}^{J}, \qquad d\mu(\v)=\sum_{j=1}^J w_j \delta(\v-\v_j), \end{equation}
with $\v_j$ denoting the chosen velocities and $w_j$ the corresponding weights. Due to this choice of $V$ the kinetic equation \eqref{eq:kin_eq_short} breaks up into a system of $J$ coupled partial differential equations,
\begin{equation} \label{eq:system_part_diff_eq} \partial_t\f^\epsi_j + \v_j\cdot\nabla_{\x}\cdot\f^\epsi_j = \frac{1}{\epsi}\left(\Maxj(\u^\epsi) - \f^\epsi_j\right), \qquad  1 \le j \le J, \end{equation}
in which $\f^\epsi_j(\x,t) \equiv \f^\epsi(\x,\v_j,t)$, and the only coupling between different velocities is through the computation of $\u^\epsi$. As $\epsi \rightarrow 0$, a Chapman-Enskog expansion allows to write
\begin{equation}\label{eq:chapman-enskog} \f^\epsi_j=\Maxj(\u^\epsi) +O(\epsi), \end{equation}
so that, injecting it in  \eqref{eq:system_part_diff_eq} and taking the mean value over $V$, we get
\begin{equation*} \partial_t \langle\Maxj(\u^\epsi)\rangle+\nabla_{\x}\cdot\langle \v_j\cdot\Maxj(\u^\epsi)\rangle=O(\epsi).\end{equation*}
Finally, the compatibility conditions \eqref{eq:maxwellian_conditions} imply
\begin{equation} \label{eq:approx} \partial_t \u^\epsi+\sum_{d=1}^D \partial_{x_d}\F^d(\u^\epsi)=O(\epsi). \end{equation}

\begin{rem} \label{rem:J2}
In the relaxation problem~\eqref{eq:kin_eq}, the minimal number $J$ of discrete velocities depends on the spatial dimension of the problem. In particular, $J=2$ in 1D and $J=4$ in 2D. Since the computational complexity of the relaxation system is proportional to $J$, we will not consider higher values of $J$ here. Hence, the one-dimensional stability analysis in section~\ref{sec:stability} will be performed specifically for $J=2$, and in the experiments in section \ref{sec:results}, we will only use $J=2$ in 1D and $J=4$ in 2D. Finally, we stress that when choosing $J=2$ in 1D, the form of the discrete kinetic system given in equation \eqref{eq:system_part_diff_eq} coincides precisely with the relaxation system introduced by Jin and Xin in \cite{jin1995relaxation}.
\end{rem}

\subsection{One-dimensional examples}
In one space dimension, we write equation \eqref{eq:cons_law_diff_system} as
\begin{equation} \label{eq:cons_law_diff_1D} \partial_t\u + \partial_x\F(\u)= 0, \end{equation}
in which $t\geq0$ (resp. $x \in \R$) represents the time (resp. space) variable, ${\u}(x,t)=(u_m(x,t))_{m=1}^M\in\R^{M}$ embodies the conserved quantities, and $\F(\u)=(F_m(\u))_{m=1}^M \in~\R^M$ represents the flux functions. Correspondingly, the kinetic equation \eqref{eq:kin_eq} becomes
\begin{equation} \label{eq:kin_eq_1d} \partial_t\f^\epsi + v \partial_x\f^\epsi = \frac{1}{\epsi}\left(\Mav(\u^\epsi) - \f^\epsi\right),\end{equation}
with the particle distribution function $\f^\epsi(x,v,t) =(f_m^\epsi(x,v,t))_{m=1}^M \in \R^M$, and the particle velocities represented as $v \in V \subset \R$.

In equation~\eqref{eq:kin_eq_1d}, the Maxwellian $\Mav$ is as yet not completely defined, as we only require the conditions~\eqref{eq:maxwellian_conditions} to be satisfied.
A physically relevant Maxwellian that corresponds to these conditions, is
\begin{equation} \label{eq:maxwellian2_1D} \Mav(\u^\epsi) = \u^\epsi + v\F(\u^\epsi), \end{equation}
in which $V=\R$ with measure
\begin{equation} \label{eq:V_measure} d\mu(v) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{v^2}{2\sigma^2}\right)dv. \end{equation}
To satisfy the subcharacteristic condition given in \eqref{eq:subchar}, the variance $\sigma$ of the measure $d\mu(v)$ needs to be chosen appropriately. For instance, in the scalar case, we require $\sigma\ge|F'(u)|$. Based on \eqref{eq:V_measure}, we choose a discrete measured velocity space with an even number $J$ of velocities that satisfy $v_{J-j+1}\equiv -v_j$.
From the measure given in \eqref{eq:V_measure}, these discrete velocities $(v_j)_{j=1}^J$ and weights $(w_j)_{j=1}^J$ can be derived as the nodes and weights of the corresponding Gauss-Hermite quadrature. For $J=2$, this results in $v_{j}=\pm\sigma$, with corresponding weights $w_j=1/2$.

An alternative suggestion, that is equivalent to the above choice for $J=2$ and $\sigma=1$, was proposed in \cite{bouchut1999construction}. There, a Maxwellian of the form
\begin{equation} \label{eq:maxwellian_1D} \Mav(\u^\epsi) = \u^\epsi + \frac{\F(\u^\epsi)}{v} \end{equation}
is proposed. For both choices \eqref{eq:maxwellian2_1D} and \eqref{eq:maxwellian_1D}, the conditions \eqref{eq:maxwellian_conditions} can be seen to be satisfied. In case of the Maxwellian given in equation \eqref{eq:maxwellian_1D}, the specific values of the velocities $v_j$ need to be chosen such that the subcharacteristic condition \eqref{eq:subchar} is satisfied.  When we further restrict to a scalar case, i.e., $M=1$,
\begin{equation} \label{eq:cons_law_diff_1D_scalar} \partial_t u + \partial_{x} F(u)= 0, \end{equation}
it can be checked that the subcharacteristic condition is always satisfied as soon as
\begin{equation} \label{eq:subchar_cond} 1+\partial_u F(u)/v_j \ge 0, \quad \text{ or, equivalently, } \quad |v_j|\ge |F'(u)|, \quad 1\le j \le J. \end{equation}
Note again that all boldfaced typesetting is removed for a scalar case.

For the numerical illustrations, we choose concretely the following examples:
\begin{ex} The linear scalar advection equation,
	\begin{equation} \label{eq:linear_flux_1D_scalar} F(u)=a\cdot u, \qquad a \in \R. \end{equation}
\end{ex}
\begin{ex} The scalar Burgers' equation,
	\begin{equation} F(u)=u^2/2. \end{equation}
\end{ex}
\begin{ex} The one-dimensional Euler equations,
	\begin{eqnarray}
		\u &=& (\rho,\rho\bar{v},\E), \label{eq:euler_conserved_vars}\\
		\F(\u) &=& (\rho\bar{v},\rho\bar{v}^2+\pr, \E+\pr\bar{v}),
	\end{eqnarray}
	with the equation of state
	\begin{equation} \pr =  (\gamma-1)\left(\E-\frac{1}{2}\rho\bar{v}^2\right). \end{equation}
\end{ex}

\subsection{Two-dimensional examples} \label{sec:2D}
In two space dimensions, we write equation \eqref{eq:cons_law_diff_system} as
\begin{equation} \label{eq:cons_law_diff_2D} \partial_t\u + \partial_{x}\F^x(\u)+\partial_{y}\F^y(\u)= 0, \end{equation}
where $t \ge 0$ (resp. $x,y \in \R$) represents the time (resp. space) variables ${\u}(x,y,t)= (u_m(x,y,t))_{m=1}^M\in~\R^{M}$ represents the conserved quantities, and $\F^{x,y}(\u)=(F^{x,y}_m(\u))_{m=1}^M\in\R^M$ represents the fluxes in the $x$ and $y$ direction respectively. Correspondingly, the kinetic equation \eqref{eq:kin_eq} becomes
\begin{equation} \label{eq:kin_eq_multidim} \partial_t\f^\epsi + v^x \partial_x \f^\epsi + v^y \partial_y \f^\epsi = \frac{1}{\epsi}(\Maxv(\u^\epsi) - \f^\epsi), \end{equation}
in which the particle distribution function $\f^\epsi(x,y,v^x,v^y,t)= (f_m^\epsi(x,y,v^x,v^y,t))_{m=1}^M \in \R^M$ and the particle velocities $\v=(v^x,v^y) \in V \subset \R^2$, with $v^{x,y}$ the velocity of the particles in the $x$ and $y$ direction respectively.

Compared to the one-dimensional setting, the choice of the Maxwellian and the description of the discrete velocity space are considerably more elaborate, and many options have been documented, see, e.g., \cite{aregba2000discrete,mieussens,bobylev}. In the numerical examples in this paper, we choose the orthogonal velocities method, see, e.g., \cite{aregba2000discrete}, which we now detail for the scalar case ($M=1$). In this method, we choose a set of velocities with varying length and direction. Specifically, we first fix a maximal velocity length $v_\textrm{max}$. We then consider $R$ different velocity lengths
\[ \rho_r = \dfrac{r}{R}v_\textrm{max}, \qquad 1\le r \le R, \]
and $4S$ different velocity directions
\[ \theta_s = \frac{s}{S}\frac{\pi}{2}, \qquad 1\le s \le 4 S, \]
with $R,S\ge 1$.
We then obtain $J=4R\,S$ velocities $\v_j=(v^x_j,v^y_j)$, $1\le j \le J$, by assigning an index $j=\left(r-1\right)4S+s$ to every length-direction pair $(r,s)$, and writing
\begin{equation}\label{eq:2D_velocity_distr} v^x_j = \rho_r\cos\left(\theta_s\right), \qquad v^y_j = \rho_r\sin\left(\theta_s\right), \qquad 1 \le r \le R, \quad 1 \le s \le 4S. \end{equation}

The Maxwellian function $\Mj$ for the $j^{th}$ equation of system \eqref{eq:system_part_diff_eq} is then chosen as:
\begin{equation} \label{eq:2D_maxwellian} \Mj(u) = u +  v^x_j\frac{F^x(u)}{a_x^2} + v^y_j\frac{F^y(u)}{a_y^2}, \end{equation}
with scalars $a_{x,y}^2 = \left\langle\left(v^{x,y}\right)^2\right\rangle$. It can be shown that for the orthogonal velocities method we have $a_x^2 = a_y^2$.

The generalization to $M>1$ is straightforward. In \cite{aregba2000discrete} it is proven that, for stability reasons, one should choose $v_{\textrm{max}}$ as follows:
\begin{equation} \label{eq:2D_v_max} v_{\textrm{max}}^2 \ge \frac{12R^2\left(\|\partial_\u \F^x\|^2 + \|\partial_\u \F^y\|^2\right)}{(R+1)(2R+1)}, \end{equation}
where $\|\cdot\|$ is the matrix norm associated with the classical 2-norm when $M > 1$.

\begin{rem} \label{rem:RS1}
As pointed out in Remark \ref{rem:J2}, in the experiments in section \ref{sec:results} we will always choose $J=4$ velocities in 2D. This is accomplished by setting $R=S=1$ in the orthogonal velocities method.
\end{rem}

\section{Projective integration} \label{sec:proj_int}
In this section, we construct a fully explicit, arbitrary order time integration method for the stiff system \eqref{eq:system_part_diff_eq}. The asymptotic-preserving property \cite{jin1999efficient} implies that, in the limit when $\epsi$ tends to zero, an $\epsi$-independent time step constraint, of the form $\Delta t = O(\Delta x)$, can be used, as the hyperbolic CFL constraint for the limiting equation \eqref{eq:cons_law_limit}. To achieve this, we will use a projective integration method \cite{gear2003projective,lafitte2012asymptotic}, which combines a few small time steps with a naive (\emph{inner}) timestepping method, such as a direct forward Euler discretization, with a much larger (\emph{projective, outer}) time step. The idea is sketched in figure \ref{fig:proj_int}.

\begin{figure}
	\begin{center}
		\figname{PI_sketch}
		\begin{tikzpicture}[node distance = 2cm, auto]

		\draw [axis] (0,0) -- (9,0);
		\draw [axis] (0,0) -- (0,3);
		\node[align=center, right] at (9,0) (time) {time};

		\draw [thick] (0.3,0.3) to [out=45,in=145] (8,0.5);

		\draw[red_dot] (0.5,2) circle (.5ex);
		\draw [dotted] (0.5,0) -- (0.5,2); \node[below] at (0.5,0) (time) {$t^{n-1}$};
		\draw[red_dot] (0.8,1.25) circle (.5ex); \draw[red_dot] (1.1,0.98) circle (.5ex); \draw[red_dot] (1.4,1.15) circle (.5ex);
		\draw [dashed] (1.1,0.98) -- (4,2.8);

		\draw[red_dot] (4,2.8) circle (.5ex);
		\draw [dotted] (4,0) -- (4,2.8); \node[below] at (4,0) (time) {$t^{n}$};
		\draw[red_dot] (4.3,2.1) circle (.5ex); \draw[red_dot] (4.6,1.83) circle (.5ex); \draw[red_dot] (4.9,1.77) circle (.5ex);
		\draw [dashed] (4.6,1.83) -- (7.5,1.4);

		\draw[red_dot] (7.5,1.4) circle (.5ex);
		\draw [dotted] (7.5,0) -- (7.5,1.4); \node[below] at (7.5,0) (time) {$t^{n+1}$};

		\path
		([shift={(-5\pgflinewidth,-5\pgflinewidth)}]current bounding box.south west)
		([shift={( 2\pgflinewidth, 2\pgflinewidth)}]current bounding box.north east);
		\end{tikzpicture}
	\end{center}
  	\caption{\label{fig:proj_int} Sketch of the general idea of a projective integration method. At each time instance ($t^n$), an explicit method is applied over a number of small time steps (red dots) so as to stably integrate the fast modes. As soon as these modes are sufficiently damped the solution is extrapolated using a much larger time step (dashed lines). }
\end{figure}

The inner and outer integrator can be selected independently. In section \ref{sec:inner}, we discuss the inner integrator. Afterwards, in section \ref{sec:outer}, we consider the outer integrator, before studying their numerical properties in section \ref{sec:stability}.

\subsection{Inner integrators}\label{sec:inner}
We intend to integrate \eqref{eq:system_part_diff_eq} on a uniform, constant in time, periodic spatial mesh with spacing $\Delta x$, consisting of $I$ mesh points $x_i=i\Delta x$, $1\le i \le I$, with $I\Delta x=1$, and a uniform time mesh with time step $\delta t$, i.e., $t^k=k\delta t$. (Generalization to space-time adaptive grids is of course straightforward.) The numerical solution on this mesh is denoted as $\f_{i,j}^k$, where we have dropped the dependence on $\epsi$ for conciseness. After discretizing in space, we obtain a semi-discrete system of ordinary differential equations
\begin{equation}\label{eq:semidiscrete} \dot{\f} = \Dt(\f),  \qquad \Dt(\f):=- \Dxv(\f) + \frac{1}{\epsi}\left(\Maxv(\u) - \f\right), \end{equation}
where $\Dxv(\cdot)$ represents a suitable discretization of the first order spatial derivative $v\partial_x$ (e.g., upwind differences).

As inner integrator, we choose an explicit scheme, for which we will, later on, use the shorthand notation
\begin{equation} \label{scheme_timestepper} \f^{k+1} =S_{\delta t}(\f^{k}),\qquad k = 0, 1, \ldots \end{equation}
The forward Euler (FE) method and Runge-Kutta methods immediately come to mind.

\paragraph{Forward Euler (FE).} The simplest time discretization routine is the forward Euler method,
\begin{equation} \label{eq:FE_scheme} \f^{k+1} = \f^k + \delta t\;\Dt(\f^k). \end{equation}

\paragraph{Higher-order Runge-Kutta methods.} \label{sec:inner_RK2}
To obtain higher-order accuracy in time in the inner integrator, one could also employ any Runge-Kutta method \cite{wanner1991solving,Hairer1993}, such as the second order method
\begin{align} \label{eq:RK2_scheme}
	\mathbf{k}_1 &= \Dt(\f^k), \\
  	\mathbf{k}_2 &= \Dt\left(\f^k + \frac{\delta t}{2}\mathbf{k}_1\right), \\
	\f^{k+1}& = \f^k + \delta t\;\mathbf{k}_2.
\end{align}
How to choose the inner integrator will be discussed in section \ref{sec:stability}.

\subsection{Outer integrators}\label{sec:outer}
In equation \eqref{eq:system_part_diff_eq}, the small parameter $\epsi$ in the relaxation term leads to the classical time step restriction of the form $\delta t = O(\epsi)$ for the inner integrator. However, as $\epsi$ goes to $0$, we obtain the limiting equation \eqref{eq:approx} for which a standard finite volume/forward Euler method only needs to satisfy a stability restriction of the form $\Delta t \le C\Delta x$, with $C$ a constant that depends on the specific choice of the scheme and the parameters of the equation.

In \cite{lafitte2012asymptotic}, it was proposed to use a projective integration method to accelerate such a brute-force integration; the idea, originating from \cite{gear2003projective}, is the following. Starting from a computed numerical solution $\f^n$ at time $t^n=n\Delta t$, one first takes $K+1$ \emph{inner} steps of size $\delta t$,
\begin{equation} \f^{n,k+1} =S_{\delta t}(\f^{n,k}),\qquad 0 \le k \le K, \end{equation}
in which the superscript pair $(n,k)$ represents the numerical solution at $t^{n,k}=n\Delta t +k\delta t$. The aim is to obtain a discrete derivative to be used in the \emph{outer} step to compute $\f^{n+1}=\f^{n+1,0}$ via extrapolation in time, e.g.,
\begin{equation} \label{PFE_scheme} \f^{n+1} = \f^{n,K+1} + \left(\Delta t - (K + 1)\delta t\right)\frac{\f^{n,K+1} - \f^{n,K}}{\delta t}. \end{equation}
This method is called projective forward Euler, and it is the simplest instantiation of this class of integration methods \cite{gear2003projective}.

In \cite{lejon13}, a higher-order projective integration method was constructed by replacing each time derivative evaluation $\mathbf{k}_s$ in a classical Runge-Kutta method by $K+1$ steps of an inner integrator as follows (with $\f^{n,0}=\f^n$ for consistency):
\begin{align}
	s = 1 : & \begin{cases}
		\f^{n,k+1} &= \f^{n,k} + \delta t\;\Dt(\f^{n,k}), \qquad 0 \le k \le K \\
		\mathbf{k}_1 &= \dfrac{\f^{n,K+1}-\f^{n,K}}{\delta t}
	\end{cases}\\
	2 \le s \le S :& \begin{cases}
		\f^{n+c_s,0}_s&=\f^{n,K+1} + (c_s\Delta t-(K+1)\delta t) \sum_{l=1}^{s-1}\dfrac{a_{s,l}}{c_s} \mathbf{k}_l, \\
		\f^{n+c_s,k+1}_s&=  \f^{n+c_s,k}_s + \delta t\;\Dt(\f^{n+c_s,k}_s), \qquad 0 \le k \le K \\
		\mathbf{k}_s &= \dfrac{\f^{n+c_s,K+1}_s-\f^{n+c_s,K}_s}{\delta t}
	\end{cases}\\
	& \f^{n+1} = \f^{n,K+1} + (\Delta t-(K+1)\delta t)\sum_{s=1}^{S}b_s \mathbf{k}_s.
\end{align}
To ensure consistency, the RK matrix $\mathbf{a}=(a_{s,l})_{s,l=1}^S$, weights $\mathbf{b}=(b_s)_{s=1}^S$, and nodes $\mathbf{c}=(c_s)_{s=1}^S$ satisfy (see, e.g., \cite{Hairer1993}) the conditions $0\le b_s \le 1$ and $0 \le c_s \le 1,$ as well as
\begin{align} \sum_{s=1}^Sb_s=1,\qquad \sum_{l=1}^{S-1} a_{s,l} =c_s, \quad 1 \le s \le S.&\label{eq:cond3} \end{align}
(Note that these assumptions imply that $c_1=0$ by the convention that $\sum_{1}^0\cdot=0$.)

In the numerical experiments, we will specifically use projective Runge-Kutta methods of orders 2 and 4, represented by the Butcher tableaux in figure \ref{table:butcher}.

\begin{figure}[t]
	\begin{center}
		\figname{butcher_tableaux}
		\begin{tikzpicture}
		\node(gen0) at(0,0){};
		\node(genx) at(0.8,0){};
		\node(geny) at(0,1){};
		\node(c) at (-0.2,0.5){$\mathbf{c}$};
		\node(A) at (0.3,0.5){ $\mathbf{a}$};
		\node(bT) at (0.3,-0.2){$\mathbf{b}^T$};
		\draw[](-0.5,0) -- (genx);
		\draw[](0,-0.5) -- (geny);
		\node(sec0) at (5,0){};
		\node(secx) at (7,0){};
		\node(secy) at (5,1){};
		\draw[] (4.45,0) -- (secx);
		\draw[] (5,-0.5) --(secy);
		\node(c1_sec) at (4.6,0.7) {$0$};
		\node (c2_sec) at (4.6,0.2){$1/2$};
		\node (a21_sec) at (5.4,0.2) {$1/2$};
		\node (b1_sec) at (5.4,-0.2) {$0$};
		\node (b2_sec) at (6.4, -0.2) {$1$};
		\node(fourth0) at (10,0){};
		\node (fourthx) at (14,0){};
		\node (fourthy) at (10,2){};
		\draw[] (9.45,0) -- (fourthx);
		\draw[] (10,-0.5) -- (fourthy);
		\node(b1_fourth) at (10.3,-0.2) {$1/6$};
		\node (b2_fourth) at (11.3,-0.2) {$1/3$};
		\node (b3_fourth) at (12.3,-0.2) {$1/3$};
		\node (b4_fourth) at (13.3,-0.2) {$1/6$};
		\node (a41_fourth) at (10.3,0.2) {$0$};
		\node (a42_fourth) at (11.3,0.2) {$0$};
		\node (a43_fourth) at (12.3,0.2) {$1$};
		\node (a31_fourth) at(10.3,0.7) {$0$};
		\node (a32_fourth) at (11.3,0.7) {$1/2$};
		\node (a21_fourth) at (10.3,1.2) {$1/2$};
		\node (c4_fourth) at (9.6,0.2) {$1$};
		\node (c3_fourth) at (9.6,0.7) {$1/2$};
		\node (c2_fourth) at (9.6,1.2) {$1/2$};
		\node (c1_fourth) at (9.6,1.7) {$0$};
		\end{tikzpicture}
	\end{center}
	\caption {\label{table:butcher} Butcher tableaux for Runge-Kutta methods. Left: general notation; middle: RK2 method (second order); right: RK4 method (fourth order).}
\end{figure}

\subsection{Stability of projective integration}\label{sec:stab_pi}
We now briefly discuss the main stability properties of projective Runge-Kutta methods as derived in \cite{lejon13}. To this end, we introduce the test equation and its corresponding inner integrator,
\begin{equation}\label{eq:test_eq} \dot{y}=\lambda y,  \qquad y^{k+1}=\tau(\lambda\delta t)y^k, \qquad \lambda \in \C. \end{equation}
As in \cite{gear2003projective}, we call $\tau(\lambda\delta t)$ the \emph{amplification factor} of the inner integrator. (For instance, if the inner integrator is the forward Euler scheme, we have $\tau(\lambda\delta t)=1+\lambda\delta t$.) The inner integrator is stable if $\left|\tau\right|\le 1$. The question then is for which subset of these values the projective integration method is also stable.

Considering projective forward Euler, it can easily be seen from \eqref{PFE_scheme} that the projective forward Euler method is stable if
\begin{equation}\label{e:stab_cond} \left|\left[\left(\dfrac{\Delta t-(K+1)\delta t}{\delta t}+1\right)\tau - \dfrac{\Delta t-(K+1)\delta t}{\delta t}\right]\tau^K \right|\le 1, \end{equation}
for all eigenvalues $\tau$ of the inner integrator for the kinetic equation \eqref{eq:kin_eq_1d}. The goal is to take a projective time step $\Delta t = O(\Delta x)$, whereas $\delta t=O(\epsi)$ necessarily to ensure stability of the inner brute-force forward Euler integration.  Since we are interested in the limit $\epsi\to 0$ for fixed $\Delta x$, we look at the limiting stability regions as $\Delta t/\delta t\to \infty$.  In this regime, it is shown in
\cite{gear2003projective} that the values $\tau$ for which the condition \eqref{e:stab_cond} is satisfied lie in the union of two separated disks $\disk^{PFE}_1\cup\mathcal{D}^{PFE}_2$ where
\begin{equation}\label{eq:stab_pfe}
	\disk^{PFE}_1=\disk\left(1-\dfrac{\delta t}{\Delta t},\dfrac{\delta t}{\Delta t}\right)\text{ and }
	\disk^{PFE}_2=\disk\left(0,\left(\dfrac{\delta t}{\Delta t}\right)^{1/K}\right),
\end{equation}
and $\disk(\kappa,\mu)$ denotes the disk with center $(\kappa,0)$ and radius $\mu$. One then aims at positioning the eigenvalues that correspond to modes that are quickly damped by the time-stepper in $\disk^{PFE}_2$, whereas the eigenvalues in $\disk_1^{PFE}$ should correspond to slowly decaying modes.  The projective integration method then allows for accurate integration of the modes in $\disk^{PFE}_1$ while maintaining stability for the modes in $\disk^{PFE}_2$.

We have the following result that compares the stability regions of higher-order projective Runge-Kutta methods to those of projective forward Euler in the limit when $\delta t/\Delta t$ tends to $0$ \cite{lejon13}.
\begin{thm}[Stability of higher-order projective Runge-Kutta methods]\label{thm:stab}
	Assume the inner integrator is stable, i.e., $\left|\tau\right|\le1$, and $K$ and $\Delta t$ are chosen in such a way that the projective forward Euler 			method is stable. Then, any projective Runge-Kutta method satisfying the conditions \eqref{eq:cond3}, as well as the convexity condition
	\begin{equation}\label{eq:rk_convex} 0 \le a_{s,l} \le c_s, \qquad 1\le l \le s,\;\; 1\le s\le S, \end{equation}
	is also stable.
\end{thm}
Such a result is classical for regular Runge-Kutta methods. In \cite{lejon13}, the proof is given in the projective Runge-Kutta case, revealing that the above property holds both for the stability domain corresponding to slow eigenvalues and for the stability domain corresponding to quickly damped eigenvalues.

Additionally, it is shown that in the limit when $\delta t/\Delta t$ tends to $0$, the stability region breaks up into two regions $\mathcal{R}_{1}^{PRK}$ and $\mathcal{R}_{2}^{PRK}$ that satisfy
\[ \mathcal{R}^{PRK,q+1}_1 \supseteq \mathcal{R}^{PRK,q}_1 \supseteq \disk^{PFE}_1 \text{ and }\mathcal{R}^{PRK,q+1}_2 \supseteq \mathcal{R}^{PRK,q}_2 \supseteq \disk^{PFE}_2, \qquad \forall q \in \{1,2,...\}, \]
in which the constant $q$ indicates the order of the specific Runge-Kutta method. The proof relies on asymptotic expansions that are similar to (but much more tedious than) those mentioned in \cite{gear2003projective}, see \cite{lejon13}.

The main conclusion is that, whereas the stability regions of higher-order projective Runge-Kutta methods differ from those of projective forward Euler in their precise shape, their qualitative dependence on the parameters of projective integration ($\delta t$, $K$ and $\Delta t$) is identical, and method parameters that are suitable for projective forward Euler will also be suitable for the higher-order projective Runge-Kutta method.

\section{Numerical properties} \label{sec:stability}
The parameters $\delta t$, $K$, and $\Delta t$ of the projective integration method can be determined by imposing that all the eigenvalues of the selected inner integrator scheme fall into the stability region of the projective integration method. Therefore, in section \ref{subsec:inner_integrator_spectrum}, we will calculate the spectrum of the inner integrators constructed for the relaxation system given in \eqref{eq:semidiscrete}, and this in the specific case $J=2$ (see Remark \ref{rem:J2}). The general result for any (even) number of velocities $J$ can be found in \cite{lejon13}. For $J=2$, we are able to derive explicit asymptotic expansions for both the fast and slow eigenvalues in the spectrum. Afterwards, we derive suitable choices for the projective integration parameters (section \ref{sec:param}). While the numerical experiments also deal with systems of nonlinear hyperbolic conservation laws in multiple space dimensions, the analysis is restricted to a one-dimensional, scalar, linear setting.

\subsection{Spectrum of inner integrators} \label{subsec:inner_integrator_spectrum}
To compute bounds on the spectrum of the inner integrator of the kinetic equation \eqref{eq:kin_eq} with a linear Maxwellian
\[ \Mv(u)=u+\dfrac{u}{v}, \]
we first rewrite the semi-discretized kinetic equation \eqref{eq:semidiscrete} in the (spatial) Fourier domain,
\begin{equation}\label{eq:fourier_semidiscr}
\partial_t \Fhat(\zeta) = \Bcal\;\Fhat(\zeta), \qquad
\Bcal = \dfrac{1}{\epsi}\left(-\epsi\D + \M\P - \Id\right),
\end{equation}
with $\Fhat \in \R^2$, the matrices $\Bcal$, $\M$, $\P \in \R^{2\times 2}$, and  $\Id$ the identity matrix of dimension $J=2$. In \eqref{eq:fourier_semidiscr}, the matrix $\D$ represents the (diagonal) Fourier matrix of the spatial discretization chosen for the convection part, $\P$ is the Fourier matrix of the averaging of $f$ over the positive and negative velocities,
\begin{equation} \label{eq:P_spectrum} \P:= \frac{1}{2}\left( \begin{array}{ccc} 1 & 1 \\ 1 & 1 \end{array} \right), \end{equation}
and the matrix $\M$ represents the Fourier transform of the Maxwellian,
\begin{equation} \label{eq:M_spectrum} \M = \Id +\mathbf{V}^{-1}, \end{equation}
with $\mathbf{V} = \diag([\vstar,-\vstar])$ and $\vstar$ the chosen (discrete) velocity component.

Since we are utilizing a symmetric velocity space, we have the following property on the diagonal elements of the matrix $\D$: $D_1 = \bar{D_2}$. Hence, from now on, we write the diagonal elements of $\D$ as
\[ D_{1,2} = \alpha \pm \imath\beta, \]
in which $\alpha$ and $\beta$ depend on the spatial discretization, the velocity value $\vstar$ and the Fourier mode $\zeta$.

Since in the Fourier domain we are calculating the eigenvalues of $(2\times 2)$ matrices, it is easy to prove the following theorem, which is a specific case of the spectrum obtained in \cite{lejon13}.
\begin{thm} \label{thm:spectrum}
	Under the above assumptions, the spectrum of the matrix $\Bcal=\dfrac{1}{\epsi}\left(\M\P - \Id - \epsi\D\right)$ contains one slow eigenvalue $\lambda_1$ and one fast eigenvalue $\lambda_2$ which can be expanded as,
	\begin{equation} \label{eq:spectrum_B}
		\begin{aligned}
		\lambda_1 &= -\alpha + \beta^2\left(\frac{1 - \vstar^2}{\vstar^2}\right)\epsi + O(\epsi^3) + \imath\left(-\frac{\beta}{\vstar} + 2\frac{\beta^3}{\vstar}\left(\frac{1-\vstar^2}{\vstar^2}\right)\epsi^2 + O(\epsi^4)\right)& \\
		\lambda_2 &= -\frac{1}{\epsi} - \alpha - \beta^2\left(\frac{1 - \vstar^2}{\vstar^2}\right)\epsi + O(\epsi^3) + \imath\left(\frac{\beta}{\vstar} - 2\frac{\beta^3}{\vstar}\left(\frac{1-\vstar^2}{\vstar^2}\right)\epsi^2 + O(\epsi^4)\right).&
		\end{aligned}
	\end{equation}
	Consequently, the spectrum of $\Bcal$ can be written as,
	\begin{equation}
	\mbox{\emph{Sp($\Bcal$)}}
	\subset
	\disk\left(-\dfrac{1}{\epsi}, \max_{\zeta}\left(\sqrt{\alpha^2+\frac{\beta^2}{\vstar^2}}\right)\right)
	\cup
	\{\lambda_1\}.
	\end{equation}
\end{thm}
\begin{proof}
Since we discretized in space, the Fourier mode $\zeta$ is a discrete variable, $\zeta_i=i2\pi\Delta x$, $i=1, ..., I$. The matrix product $\M\P$ and matrix $\D(\zeta)$ in the Fourier domain become:
\begin{equation} \label{eq:MP_and_D}
\M\P = \frac{1}{2}\left( \begin{array}{ccc} 1+1/\vstar & 1+1/\vstar \\ 1-1/\vstar & 1-1/\vstar \end{array} \right), \qquad
\D = \left( \begin{array}{ccc} z^{+} & 0 \\ 0 & z^{-} \end{array} \right),
\end{equation}
where we used the matrices $\M$ and $\P$ given in equations \eqref{eq:M_spectrum} and \eqref{eq:P_spectrum}. Furthermore, in \eqref{eq:MP_and_D} we adopted the shorthand notation $z^{\pm} = \alpha \pm \imath\beta$, where $\alpha$ and $\beta$ depend on the spatial discretization method, the velocity value $\vstar$ and the Fourier mode $\zeta$.

We start by calculating the spectrum of the matrix $\Acal = \epsi\Bcal$. The eigenvalues $\tilde\lambda$ of $\Acal$ can be obtained as the roots of its characteristic polynomial $\chi_{\Acal}(\tilde\lambda)$, which is given by,
\begin{equation} \label{eq:Acal_char_poly}
	\begin{aligned}
	\chi_{\Acal}(\tilde\lambda) &= |\Acal - \tilde\lambda\Id|& \\
	&= \tilde\lambda^2 + (1 + \epsi(z^{+}+z^{-}))\tilde\lambda + \epsi\left(\frac{1+\vstar}{2\vstar}z^{+}-\frac{1-\vstar}{2\vstar}z^{-}\right) + \epsi^2z^{+}z^{-} = 0.&
	\end{aligned}
\end{equation}
Using the equalities,
\[ z^{+}+z^{-} = 2\alpha, \qquad z^{+}-z^{-} = \imath2\beta, \qquad z^{+}z^{-} = \alpha^2+\beta^2, \]
the roots of equation \eqref{eq:Acal_char_poly} can be calculated as:
\begin{equation} \label{eq:Acal_eigs} \tilde\lambda_{1,2} = \frac{-1-2\alpha\epsi \pm \sqrt{D}}{2}, \end{equation}
where the discriminant $D$ of equation \eqref{eq:Acal_char_poly} can be written as,
\[ D = 1 - 4\beta^2\epsi^2 - \imath4\frac{\beta}{\vstar}\epsi. \]
Employing a Taylor series expansion for $\sqrt{D}$ for $\epsi \rightarrow 0$ in equation \eqref{eq:Acal_eigs} leads to:
\begin{equation*} \sqrt{D} = 1 + 2\beta^2\left(\frac{1 - \vstar^2}{\vstar^2}\right)\epsi^2 + O(\epsi^4) + \imath\left(-2\frac{\beta}{\vstar}\epsi + 4\frac{\beta^3}{\vstar}\left(\frac{1-\vstar^2}{\vstar^2}\right)\epsi^3 + O(\epsi^5)\right). \end{equation*}
Plugging this expansion into equation \eqref{eq:Acal_eigs} we obtain,
\begin{equation*}
	\begin{aligned}
	\tilde\lambda_1 &= \frac{-1-2\alpha\epsi + \sqrt{D}}{2}& \\
	& = -\alpha\epsi + \beta^2\left(\frac{1 - \vstar^2}{\vstar^2}\right)\epsi^2 + O(\epsi^4) + \imath\left(-\frac{\beta}{\vstar}\epsi + 2\frac{\beta^3}{\vstar}\left(\frac{1-\vstar^2}{\vstar^2}\right)\epsi^3 + O(\epsi^5)\right)& \\
	\tilde\lambda_2 &= \frac{-1-2\alpha\epsi - \sqrt{D}}{2}& \\
	& = -1 - \alpha\epsi - \beta^2\left(\frac{1 - \vstar^2}{\vstar^2}\right)\epsi^2 + O(\epsi^4) + \imath\left(\frac{\beta}{\vstar}\epsi - 2\frac{\beta^3}{\vstar}\left(\frac{1-\vstar^2}{\vstar^2}\right)\epsi^3 + O(\epsi^5)\right).&
	\end{aligned}
\end{equation*}
Finally, the dominant and fast eigenvalues, $\lambda_1$ and $\lambda_2$, of the matrix $\Bcal = \dfrac{1}{\epsi}\Acal$ are then given by
\begin{equation*}
	\begin{aligned}
	\lambda_1 &= -\alpha + \beta^2\left(\frac{1 - \vstar^2}{\vstar^2}\right)\epsi + O(\epsi^3) + \imath\left(-\frac{\beta}{\vstar} + 2\frac{\beta^3}{\vstar}\left(\frac{1-\vstar^2}{\vstar^2}\right)\epsi^2 + O(\epsi^4)\right)& \\
	\lambda_2 &= -\frac{1}{\epsi} - \alpha - \beta^2\left(\frac{1 - \vstar^2}{\vstar^2}\right)\epsi + O(\epsi^3) + \imath\left(\frac{\beta}{\vstar} - 2\frac{\beta^3}{\vstar}\left(\frac{1-\vstar^2}{\vstar^2}\right)\epsi^2 + O(\epsi^4)\right).&
	\end{aligned}
\end{equation*}
\end{proof}

When we write the Fourier transform of the inner forward Euler scheme \eqref{eq:FE_scheme} as
 \begin{equation}\label{eq:fourier_fe} \Fhat^{k+1} = \mathbf{S}_{\delta t}\;\Fhat^{k} = \left(\Id+\delta t\Bcal \right)\Fhat^k, \end{equation}
it is clear that the amplifications factors $\boldsymbol{\tau}=(\tau_1,\tau_2)$ of the forward Euler scheme, which are the eigenvalues of $\mathbf{S}_{\delta t}$, and the eigenvalues $\boldsymbol{\lambda}=(\lambda_1,\lambda_2)$ of the matrix $\Bcal$ are related via
\begin{equation}\label{eq:tau_fe} \tau_j = 1 + \delta t \lambda_j, \qquad j \in \{1,2\}. \end{equation}
Thus, the spectrum of an inner forward Euler time-stepper satisfies
\begin{equation}\label{eq:spec_inner_fe}
\textrm{Sp}(\Id+\delta t\Bcal)
\subset
\disk\left(1-\dfrac{\delta t}{\epsi}, \delta t\max_{\zeta}\left(\sqrt{\alpha^2+\frac{\beta^2}{\vstar^2}}\right)\right)
\cup
\{1+\lambda_1\delta t\},
\end{equation}
with $\lambda_1$ given in Theorem \ref{thm:spectrum}.

For higher-order Runge-Kutta inner integrators (of order $Q$), we have
\begin{equation} \label{eq:tau_rk} \tau_j = 1+\sum_{q=1}^Q\frac{(\delta t \lambda_j)^q}{q!}, \qquad j \in \{1,2\}, \end{equation}
and thus
\begin{equation}\label{eq:spec_inner_rk}
\textrm{Sp}\left(\Id+\sum_{q=1}^Q\frac{(\delta t \Bcal)^q}{q!}\right)
\subset
\disk\left(1+\sum_{q=1}^Q\dfrac{(-1)^q}{q!}\left(\dfrac{\delta t}{\epsi}\right)^q, C\delta t\max_{\zeta}\left(\sqrt{\alpha^2+\frac{\beta^2}{\vstar^2}}\right)\right)
\cup
\left\{1+\sum_{q=1}^Q\frac{(\delta t \lambda_1)^q}{q!}\right\}, \end{equation}
in which the constant $C$ depends on the order of the Runge-Kutta method. The spectrum in formula \eqref{eq:spec_inner_rk} can be obtained by transforming the spectrum of $\Bcal$ in \eqref{eq:spectrum_B} by the particular expression of the Runge-Kutta inner integrator amplification factor given by \eqref{eq:tau_rk}.

\subsection{Method parameters} \label{sec:param}
In this section the projective integration method parameters will be determined by ensuring that the spectrum of the inner integrator falls within the stability region of the outer integrator. First, we select a suitable inner time step $\delta t$ such that the fast modes are quickly damped (section~\ref{sec:choose-dt}). Then, we choose the outer time step $\Delta t$ commensurate with the slow part of the evolution (section~\ref{sec:choose-Dt}). Finally, we fix $K$ to ensure overall stability (section~\ref{sec:choose-K}).

\subsubsection{Choice of inner integrator and time step}\label{sec:choose-dt}
Let us first discuss the effect of the choice of the inner integrator.  To this end, we look at the discretization error and the desired stability properties.

Concerning stability, we deduce from the stability properties of the projective integration method (see section \ref{sec:stab_pi}) that it is preferable to center the part of the spectrum of the inner time-stepper corresponding to quickly damped modes around $0$.  Since, for forward Euler, these fast modes are given by \eqref{eq:spec_inner_fe}, we choose, for inner forward Euler, $\delta t=\epsi$.

For higher-order inner integrators with even order, we immediately see that one cannot center eigenvalues corresponding to the quickly damped modes around $0$.  For instance, for a second order Runge-Kutta method we have that
\begin{equation} \label{eq:max_eig_dist} \min_{\delta t\lambda_2} \tau_2 = \min_{\delta t\lambda_2}\left(1+\delta t\lambda_2+\dfrac{1}{2}(\delta t\lambda_2)^2\right) = \dfrac{1}{2}, \end{equation}
and this value is reached for $\delta t=\epsi$. Moreover, as we will show below, the discretization error of the projective integration scheme is dominated by the error of the outer integrator, whereas the discretization error due to the inner integrator is negligible. As a consequence, we conclude that there is no point in using a higher-order time discretization for the inner integrator.

\subsubsection{Outer time step}\label{sec:choose-Dt}
Given the inner time step $\delta t=\epsi$, we can choose $\Delta t$ such that the dominant eigenvalue $\tau_1$ of the inner integrator lies inside the stability region $\mathcal{R}_1^{PRK,q}$. Let us first look at the projective forward Euler method, with stability regions \eqref{eq:stab_pfe}. We have the following condition on $\Delta t$ such that $\tau_{1}$ is contained within $\disk_1^{PFE}$,

\[ \sqrt{\left(\Re(\tau_1) -  \left(1-\frac{\epsi}{\Delta t}\right)\right)^2 + \Big(\Im(\tau_1)\Big)^2} \le \left(\frac{\epsi}{\Delta t}\right)^2. \]
Using \eqref{eq:tau_fe} and $\lambda_1$ given in theorem \ref{thm:spectrum}, we deduce a bound of the form
\[ \Delta t \le \min_\zeta\left(\frac{2\alpha \vstar^2}{\alpha^2\vstar^2 + \beta^2}\right). \]
For the projective forward Euler method and first order upwind we obtain the CFL-like bound
\[ \Delta t \le \frac{\Delta x}{\vstar}. \]

We remark that this is exactly the stability condition for the forward Euler method applied to the original hyperbolic conservation law~\eqref{eq:cons_law_diff_1D_scalar}. When we combine a higher-order upwind method with forward Euler time-stepping, the time step restriction becomes much more severe, in exactly the same way as would be the case with a direct higher-order upwind/forward Euler discretization of equation~\eqref{eq:cons_law_diff_1D_scalar}. We find, for instance,
\begin{equation*}
	\begin{aligned}
	\Delta t &\le \frac{2\pi^2\vstar(3-10\pi^2\Delta x^2)}{12}\Delta x^3, \qquad &\text{(second order upwind)}& \\
	\Delta t &\le \frac{2\pi^2\vstar(3-2\pi^2\Delta x^2)}{9}\Delta x^3, \qquad &\text{(third order upwind)}.&
	\end{aligned}
\end{equation*}
This is because, for higher-order upwind methods, the dominant eigenvalues do not lie on a circle anymore. In fact, they belong to a region that is much steeper close to 1. Since the dominant stability region of projective forward Euler is always a circle, we need to choose its radius sufficiently large such that even the steepest eigenvalues fall into this circular stability region.
When we use a higher-order projective Runge-Kutta method, we should determine the bound on $\Delta t$ using its proper dominant stability region (see, e.g., \cite{lejon13}). For instance, in case of the second order projective Runge-Kutta method, we obtain the following CFL-like bound,
\[ \Delta t \le \nu\Delta x, \]
where $\nu$ is the CFL-number which depends on the coefficients $\alpha$ and $\beta$. Notice that in both cases $\Delta t$ never depends on $\epsi$.

\subsubsection{The number $K$ of inner time steps}\label{sec:choose-K}
The only parameter that remains to be chosen is the number $K$ of inner steps such that all fast eigenvalues $\tau_2$ (see equation \eqref{eq:tau_fe} or \eqref{eq:tau_rk}) corresponding to quickly damped modes are contained within the stability region $\disk_2^{PFE}$ (see equation \eqref{eq:stab_pfe}).  Given $\delta t$, $\Delta t$ and the radius of the fast eigenvalues zone introduced in \eqref{eq:spec_inner_fe} or \eqref{eq:spec_inner_rk}, we obtain a condition for $K$ of the form,
\[ c\epsi \le \left(\dfrac{\epsi}{\Delta t}\right)^{1/K},  \qquad c = \max_{\zeta}\left(\sqrt{\alpha^2+\frac{\beta^2}{\vstar^2}}\right). \]
From this we can obtain the bound
\[ K \ge \frac{1}{1+\dfrac{\log(c)}{\log(\epsi)}} + \frac{\log(\Delta t)}{\log\left(\dfrac{1}{\epsi c}\right)}. \]
Therefore, for $\epsi$ small a safe choice is taking $K\ge 2$, uniformly in $\epsi$. The computational cost of the method is then independent of $\epsi$.

Note that, when choosing an inner Runge-Kutta method of even order, the fast eigenvalues $\tau_2$ cannot be centered around $0$, see, for instance, \eqref{eq:max_eig_dist}. In that case, we need to impose that the radius of the stability region $\disk_2^{PFE}$ (see equation \eqref{eq:stab_pfe}) is sufficiently large to encompass the fast eigenvalues zone of the inner Runge-Kutta method, yielding
\[ \tilde c \le \left(\dfrac{\delta t}{\Delta t}\right)^{1/K}, \qquad \tilde c = 1+\sum_{q=1}^Q\dfrac{(-1)^q}{q!}\left(\dfrac{\delta t}{\epsi}\right)^q + C\delta t\max_{\zeta}\left(\sqrt{\alpha^2+\frac{\beta^2}{\vstar^2}}\right). \]
Taking the logarithm of both sides and rearranging terms ultimately leads to
\[ K \ge \frac{\log\left(\dfrac{1}{\delta t}\right)}{\log\left(\dfrac{1}{\tilde c}\right)} + \frac{\log(\Delta t)}{\log\left(\dfrac{1}{\tilde c}\right)}, \]
which results in a condition of the form $K\ge\log(1/\epsi)$ since $\delta t = O(\epsi)$. Hence, using an inner integrator of even order destroys the asymptotic-preserving nature of the projective integration method.

\subsection{Consistency result}
Since the goal of the simulation is to obtain $\u(\x,t)=\langle \f(\x,\v,t)\rangle$, we define the space-time discretization error at time $t=t^N$ and grid location $\x=\x_i$ as
\begin{equation}\label{eq:error} E^N_i = \left\| \u_i^N-\u(\x_i,t^N)\right\|, \end{equation}
with $\|\cdot\|$ a vector norm of choice, and $E^N=\|\boldsymbol{E}^N\|$, with $\boldsymbol{E}^N=(E_i^N)_{i=1}^I$.

From \cite{lejon13}, we cite the following consistency result:
\begin{thm}[Consistency]
	When using an inner forward Euler time-stepper \eqref{eq:FE_scheme} with $\delta t=\epsi$ and a space discretization of order $p$, a projective 			Runge-Kutta method of order $q$ has the following discretization error:
	\begin{equation} \label{eq:local_trunc_error} E^N = O\left(\Delta t^q\right) + \epsi O\left(\Delta x^p\right) + O\left((K+1)\epsi\right) + \epsi O\left(\frac{\epsi\Delta x^p + \epsi^2)}{\Delta t}\right).\end{equation}
\end{thm}

The first term in \eqref{eq:local_trunc_error} is due to the time discretization error made in the outer Runge-Kutta integrator, whereas the next two terms are the space and time discretization error respectively due to the inner integrator. The last term is due to the estimate of the time derivative operator. We remark that, in the limit of $\Delta t$ going to $0$, for fixed $\epsi$, this last term would result in divergence. However, since the goal is to create an asymptotic-preserving scheme, valid for fixed $\Delta t$ (independent of $\epsi$), while $\epsi$ tends to $0$, this is not an issue: the last term then becomes of $O(\epsi)$, and hence negligible.

\section{Applications} \label{sec:results}
Let us now illustrate the relaxation method with projective integration on a number of example systems. We first examine the one-dimensional case. In section \ref{subsec:order}, we consider the linear advection equation, and demonstrate the spatial and temporal order of the methods.  Subsequently, we investigate nonlinear conservation laws, Burgers' equation in section \ref{subsec:burgers} and the Euler equations (Sod's shock test) in section \ref{subsec:sod}. Afterwards, we consider linear advection, the dam-break problem and the Euler equations in the two-dimensional case (sections \ref{subsec:lin_2d}--\ref{subsec:euler_2d}).

\subsection{Linear advection in 1D} \label{subsec:order}
Let us first illustrate the order of the relaxation method with projective integration, both in space and time. To this end, we consider the linear advection equation, i.e., equation \eqref{eq:cons_law_diff_1D_scalar} with the linear flux function \eqref{eq:linear_flux_1D_scalar},
\begin{equation} \label{eq:lin_adv_1D_scalar} \partial_t u + \partial_x (a\cdot u)= 0, \end{equation}
in which the macroscopic unknown function $u(x,t)$ denotes the density of the particles. We compute the solution for $t\in[0,T]$ and $x\in[0, L]$, using $a = L = 1$. We impose periodic boundary conditions and choose a smooth initial condition:
\begin{equation} \label{order_test_BC} u(x,0) = \exp(-100(x-0.5)^2). \end{equation}
We  compute the error $E$ at time $t=T$ using \eqref{eq:error}, in which we have used the $1$-norm.  (For the linear advection equation, the exact solution is known analytically.)

For the relaxation method, we use the kinetic equation \eqref{eq:kin_eq_1d}, in which we discretize the velocity space using $J=2$ velocities. Taking into account the subcharacteristic condition \eqref{eq:subchar_cond}, the velocities are chosen as: $v_1 = -|a|$ and $v_2 = |a|$. We point out that, in this purely academic first test case, we recover the Jin-Xin relaxation system (see also Remark 2.1), and additionally, the distribution corresponding to $v_1=-|a|$ vanishes in the limit of $\epsi$ tending to zero.
As the Maxwellian, we choose \eqref{eq:maxwellian_1D}, with $F(u)=a\cdot u$. The inner integrator is a space-time discretization of equation \eqref{eq:kin_eq_1d}, in which we choose the standard upwind spatial discretizations of order $1$, $2$ and $3$ with grid spacing $\Delta x$ (that will vary throughout the experiments), combined with a forward Euler time discretization with $\delta t=\epsi$ and $\epsi=10^{-8}$. The projective integration method uses $K=2$ inner steps, and an outer time step of size $\Delta t$ (that will also vary).

\paragraph{Numerical spatial order (figures \ref{fig:order_test_PFE_PRK2} and \ref{fig:order_test_PRK4})}
To illustrate the spatial order of accuracy, we calculate the error on time $T=0.02$ and vary the grid spacing $\Delta x$ as
\begin{equation} \label{order_test_dx} \Delta x = [0.05; 0.02; 0.01; 0.005; 0.002; 0.001; 0.0005; 0.0002; 0.0001], \end{equation}
and correspondingly choose the outer time step $\Delta t$ as
\begin{equation}\label{order_test_Dt}  \Delta t = O\left(\Delta x^{p/q}\right) = C_p\Delta x^{p/q}, \end{equation}
in which the constant $C_p$ should be chosen such that the projective integration method remains stable for all choices of $\Delta x$. Figure \ref{fig:order_test_PFE_PRK2} shows the error as a function of $\Delta x$.  In the leftmost figure, time integration is done using the projective forward Euler method (PFE), for which the time order $q=1$.  The constants in \eqref{order_test_Dt} are then chosen as $C_1=0.8$, $C_2 = 20$ and $C_3 = 100$.  We clearly observe the expected spatial order. This is confirmed by fitting a least squares line through the calculated error points. The slopes of these lines correspond to the numerical order which in this case were found to be 1.07, 1.98 and 2.99. These indeed lie sufficiently close to the expected spatial order. In the middle figure, the experiment is repeated using a second order projective Runge-Kutta method (PRK2, $q=2$). For spatial orders $p \in \{1,2\}$, we choose $\Delta t = 0.5\Delta x$. In that case, the first term in \eqref{eq:local_trunc_error} will be dominant and the order in space can be observed. For $p = 3$ we put $\Delta t = C_3\Delta x^{3/2}$ and choose $C_3 = 5$. On these plots, we observe that, for the third order upwind discretization, the error curves start to level off for small values of $\Delta x$. This is due to the fact that the contribution of the spatial discretization error in \eqref{eq:local_trunc_error} becomes negligible, and the $O(\epsi)$ term becomes dominant. This term is due to the approximation of the time derivative for the projective step by a finite difference approximation of the form
\begin{equation} \label{projective_der_approx} \frac{f^{n,K+1} - f^{n,K}}{\epsi}. \end{equation}
As a consequence, the numerical error cannot be decreased further than $O\left(\sqrt{\epsi_{mach}}\right)$ with $\epsi_{mach} \approx 10^{-16}$ being the machine precision. This is the reason we choose $\epsi = 10^{-8}$. When calculating the slopes of the least squares fit we now obtain 0.99, 2.00 and 2.91 which correspond to the expected spatial orders.

Finally, we repeat the experiment using a second order Runge-Kutta method as the inner integrator (see figure \ref{fig:order_test_PFE_PRK2}, right).  We note, in agreement with the observations in section \ref{sec:param} that $K$ can now no longer be chosen independently of $\epsi$.  Here, we choose $K=21$ for $p\in\{1,2\}$ and $K=23$ for $p=3$. The numerical orders are 0.99, 2.00 and 2.88. Moreover, we again observe that the error levels off, since the remaining error is due to the finite difference approximation of the time derivative, and not due to the time discretization error of the inner integrator. Note that the error curve now levels off to a value which is about $10$ times higher than in the forward Euler case, due to a less efficient damping of the fast eigenvalues, as can be seen in \eqref{eq:spec_inner_rk}. This observation is supported by looking at equation \eqref{eq:local_trunc_error} in which the value of $K$ needs to be taken $10$ times higher for RK2 as inner integrator so as to guarantee a stable functioning of the method. From these findings we conclude that it is not useful to select a higher-order inner integrator within the projective integration framework. Therefore, in what follows we will always select FE as inner integrator.

Next, we repeat this experiment using a fourth order projective Runge-Kutta method (PRK4). We choose $\Delta t = C\Delta x$ for each $p \in \{1, 2, 3\}$ and put $C$ equal to $0.4$. The result is depicted on the left hand side plot of figure \ref{fig:order_test_PRK4}, where we again see the expected behavior (numerical orders: 0.99, 2.00 and 3.00).  To avoid unphysical oscillations associated with higher-order upwind schemes, we also performed the same experiments using an essentially non-oscillatory (ENO) spatial discretization \cite{eno_weno}, which uses an adaptive stencil thus trying to avoid stencils with large variations in the solution values. The order test of the PRK4 scheme with ENO is shown on the right hand side plot of figure \ref{fig:order_test_PRK4}. The calculated numerical orders were 0.99, 1.88 and 3.13 which are in agreement with the expected spatial orders of the ENO scheme.

\begin{figure}
	\begin{center}
		\includegraphics[scale=0.47]{ordertest_space_PFE_PRK2}
	\end{center}
  	\caption{\label{fig:order_test_PFE_PRK2} Spatial order test for PFE with FE as inner integrator (left) and PRK2 with FE (middle) and RK2 (right) as inner integrators and three different spatial orders. The error is computed using the 1-norm. On each plot, the solid lines represent the calculated error whereas the dotted line shows the expected error. }
\end{figure}

\begin{figure}
	\begin{center}
		\includegraphics[scale=0.47]{ordertest_space_PRK4_u123_ENO123}
	\end{center}
	\caption{\label{fig:order_test_PRK4} Spatial order test for PRK4 with FE as inner integrator and three different spatial orders using a) upwind differences (left plot) and b) the ENO scheme (right plot). The error is computed using the 1-norm. On each plot, the solid lines represent the calculated error whereas the dotted line shows the expected error. }
\end{figure}

\begin{figure}
	\begin{center}
		\includegraphics[scale=0.47]{ordertest_time_linadv}
	\end{center}
	\caption{\label{fig:order_test_time_lin_adv} Temporal order test for the different projective integration methods under study: PFE, PRK2 and PRK4 with FE as inner integrator and upwind differences of order 3 in space comparing results for $\epsi=10^{-5}$ (left plot) and $\epsi=10^{-8}$ (right plot). The error is computed using the 1-norm. On each plot, the solid lines represent the calculated error whereas the dotted line shows the expected error. }
\end{figure}

\paragraph{Numerical time order (figure \ref{fig:order_test_time_lin_adv})}
The temporal order of the projective integration methods is demonstrated in a slightly different manner as outlined above. Now, we fix the grid spacing $\Delta x=10^{-2}$ and vary the outer time step $\Delta t$ as
\begin{equation*}  \Delta t = [0.04; 0.02; 0.01; 0.005; 0.002; 0.001; 0.0005; 0.0002; 0.0001]. \end{equation*}
We will calculate the error on time $T=0.04$. The other simulation parameters remain the same as above. The error is now calculated by taking the 1-norm of the difference between the numerical solution and the analytical solution of the (linear) semi-discretized system \eqref{eq:semidiscrete}. By doing so we take into account the discretization error in space such that we only look at the error in time. The simulations are run for PFE, PRK2 and PRK4 with FE as inner integrator and upwind differences of order 3 in space. The results can be seen in figure \ref{fig:order_test_time_lin_adv} in which we also look at the influence of the value of $\epsi$ by choosing $\epsi=10^{-5}$ (left hand side plot) and $\epsi=10^{-8}$ (right hand side plot). It is clearly indicated on the plots that for small values of $\Delta t$ the error curves level off towards the value of the dominant term of $O(\epsi)$ in expression \eqref{eq:local_trunc_error} since the other terms in \eqref{eq:local_trunc_error} are negligible. For $\epsi=10^{-5}$ the numerical orders are 1.35 and 1.99 for PFE and PRK2 respectively. In this case there were too few meaningful points to reliably estimate the numerical order of the PRK4 method. For $\epsi=10^{-8}$ the numerical orders are 1.00, 1.99 and 3.99, corresponding to the expected time order of the different methods.

\subsection{Burgers' equation in 1D}\label{subsec:burgers}
As a second example, we consider the inviscid Burgers' equation in one spatial dimension,
\begin{equation} \label{burgers_eq} \partial_t u + \partial_x\left(\frac{u^2}{2}\right) = 0. \end{equation}
We compute the solution for $t\in[0, 1]$ and $x\in[0, L]$, using $L = 2$. We impose periodic boundary conditions and consider three different initial conditions: a Gaussian pulse $u_1(x)$, a sinc wave packet $u_2(x)$ and a sine wave $u_3(x)$ given by,
\begin{equation} \label{burgers_eq_init_cond} u_1(x) = \exp(-25(x-1)^2), \quad u_2(x) = \text{sinc}(5(x - 1)), \quad u_3(x) = \sin(\pi x). \end{equation}

For the relaxation method, we use the kinetic equation \eqref{eq:kin_eq_1d} together with the more realistic form of the Maxwellian given in \eqref{eq:maxwellian2_1D} with $F(u)=u^2/2$. (Note that, compared to the previous example with linear advection, this only requires a change in one line of the code.) We discretize the velocity space using $J=2$ velocities which are obtained as the nodes of Gauss-Hermite quadrature. The inner integrator is a space-time discretization of equation \eqref{eq:kin_eq_1d}, in which we again choose the standard upwind spatial discretizations of order $1$, $2$ and $3$ with grid spacing $\Delta x$ (that will vary throughout the experiments), combined with a forward Euler time discretization with $\delta t=\epsi$ and $\epsi=10^{-8}$.  We also consider the third order ENO scheme. The projective integration method uses $K=2$ inner steps, and an outer time step of size $\Delta t$ (that will also vary).

We first perform a numerical simulation using a third order ENO spatial discretization and a fourth order projective Runge-Kutta method (PRK4) with  $\Delta x = 10^{-2}$ and $\Delta t=5\cdot 10^{-3}$. The results are shown in figure \ref{fig:burgers_PRK4_FE_ENO3}. We clearly see that the discontinuities are nicely captured without the appearance of spurious oscillations.

\begin{figure}
	\begin{center}
		\includegraphics[scale=0.47]{burgers_PRK4_FE_ENO3_3waves}
	\end{center}
	\caption{\label{fig:burgers_PRK4_FE_ENO3} Evolution of the numerical solution of Burgers' equation obtained with PRK4 with FE and the third order ENO scheme. }
\end{figure}

Next, we also investigate the temporal order of the methods. Since the analytical solution of Burgers' equation is not available explicitly, the error is computed with respect to a reference solution that is obtained using a high-order simulation of \eqref{burgers_eq} with the PRK4 method with FE as inner integrator and upwind differences of order 3 in space using a grid spacing $\Delta x = 10^{-2}$ and a time step $\delta t = 10^{-8}$. The outer time step is chosen to be $\Delta t=10^{-6}$. We again examined the influence of the value of $\epsi$ by choosing $\epsi=10^{-5}$ (left hand side plot) and $\epsi=10^{-8}$ (right hand side plot) in figure \ref{fig:order_test_time_burgers}. It is observed on the plots that for small values of $\Delta t$ the error curves level off towards the value of the dominant term of $O(\epsi)$ in expression \eqref{eq:local_trunc_error} since the other terms in \eqref{eq:local_trunc_error} are negligible. For $\epsi=10^{-5}$ the numerical orders are 1.04 and 1.96 whereas for $\epsi=10^{-8}$ the numerical orders are 1.00, 1.99 and 3.96 which are in agreement with the expected time order of the different methods. For the temporal order test the initial solution consisted of a Gauss curve centered around the middle of the domain $x\in[0,2]$ (i.e. function $u_1(x)$ in \eqref{burgers_eq_init_cond}) and the error was calculated on $t=0.04$. The results are comparable to the linear advection case, see figure \ref{fig:order_test_time_burgers}.

\begin{figure}
	\begin{center}
		\includegraphics[scale=0.47]{ordertest_time_burgers}
	\end{center}
	\caption{\label{fig:order_test_time_burgers} Temporal order test for Burgers' equation for the different projective integration methods under study: PFE, PRK2 and PRK4 with FE as inner integrator and upwind differences of order 3 in space comparing results for $\epsi=10^{-5}$ (left plot) and $\epsi=10^{-8}$ (right plot). The error is computed using the 1-norm. On each plot the solid lines represent the calculated error whereas the dotted line shows the expected error. }
\end{figure}

\subsection{Sod's shock test in 1D}\label{subsec:sod}
Sod's shock test is an important numerical test to check how a numerical method captures shock waves \cite{sod1978survey}. The test involves the Euler equations in one spatial dimension for mass, momentum and energy, $\u=(\rho,\rho\bar{v},\E)$,
\begin{equation}\label{Euler_eq}
	\begin{cases}
		\partial_t \rho + \partial_x \left(\rho\bar v\right) = 0, \\
		\partial_t \left(\rho\bar v\right) + \partial_x \left(\rho\bar v^2 + \pr\right) = 0, \\
		\partial_t \E + \partial_x \left(\E + \pr\bar v\right) = 0,
	\end{cases}
\end{equation}
in which the pressure $\pr$ is determined via the equation of state
\begin{equation}
	\label{eq_of_state}
	\pr = (\gamma - 1)\left(\E - \frac{1}{2}\rho\bar{v}^2\right),
\end{equation}
in which the constant $\gamma$ equals $7/5$ in case of a diatomic perfect gas \cite{sod1978survey}.

As an initial condition, Sod's shock test imposes
\begin{equation} \label{Euler_eq_IC}
	\begin{aligned}
		\rho(x,0) &= \left\{ \begin{array} {l@{\qquad \quad}l} 1 & x \le L/2\\ 0.125 & x > L/2 \end{array} \right.&, \\
		\pr(x,0) &=  \left\{ \begin{array} {l@{\qquad \quad}l} 1 & x \le L/2\\ 0.1 & x > L/2 \end{array} \right.&, \\
		\bar v(x,0) &= 0&.
	\end{aligned}
\end{equation}
Additionally, we impose outflow boundary conditions. We perform the simulation over the spatial domain $x \in [0,L]$ with $L=1$ and $t\in [0,0.22]$. The particular choice of the time interval allows for a clear visualization of the three different characteristic waves (see below) and an easy comparison with the available literature.

For the relaxation method, we use the kinetic equation \eqref{eq:kin_eq_1d}, in which we discretize the velocity space with $J=2$ velocities corresponding to the nodes of Gauss-Hermite quadrature. The Maxwellian is chosen as
\begin{equation} \label{eq:sod_maxwellian}
	\Mav(\textbf{u}) = \left( \begin{array}{c} \rho + v(\rho\bar v) \\ \rho\bar v + v(\rho\bar v^2 + \pr) \\ \E + v(\E + \pr)\bar v \end{array} \right). \end{equation}
The inner integrator uses a third order spatial ENO discretization with $\Delta x = 5\cdot 10^{-3}$ and a forward Euler time discretization with $\delta t = \epsi = 10^{-8}$.  As the outer method, we choose the fourth order projective Runge-Kutta method (PRK4), using $K=2$ inner steps and an outer step of size $\Delta t = 0.5\Delta x$.

The result is illustrated in figure \ref{fig:sod_PRK4_ENO3}, where we plot density, velocity, energy and pressure at time $t=0.22$, along with the analytical solution, calculated with an exact Riemann solver, for comparison.

\begin{figure}
	\begin{center}
		\includegraphics[scale=0.50]{sod2_PRK4_FE_ENO3}
	\end{center}
	\caption{\label{fig:sod_PRK4_ENO3} Evolution of several fluid properties in Sod's shock test obtained with PRK4 and FE as inner integrator and a third order ENO finite difference discretization in space. }
\end{figure}

We clearly see the development of three characteristic waves. The first wave corresponds to a rarefaction wave propagating to the left since initially the pressure and density on the left half side of the shock tube are higher than on the right side. Secondly, a contact discontinuity is observed. This right propagating wave corresponds to the initial discontinuity of \eqref{Euler_eq_IC}. Finally, a shock wave propagating to the right has also appeared. It appears that only pressure and velocity remain continuous over the contact discontinuity and exhibit a very flat state. Furthermore, it can be observed that all the quantities are discontinuous over the shock wave. The projective integration method captures all these phenomena, without developing undesired oscillations and without smoothing out the discontinuities too much.

\subsection{Linear advection in 2D}\label{subsec:lin_2d}
Let us now turn to problems in two spatial dimensions. We again start with the linear advection equation, which, in dimension 2, reads
\begin{equation} \label{2D_lin_adv} \partial_t u + a\partial_x u + b\partial_y u = 0, \end{equation}
in which $a,b \in \R$ are the constant advection speeds along the $x$- and $y$-direction respectively. The macroscopic unknown function $u(x,y,t)$ denotes the two-dimensional density of particles. In the simulations we integrate over $t\in[0,1]$ and $(x,y)\in[0,L]^2$, and set $a = b = L =1$. We impose periodic boundary conditions and start from a Gaussian pulse centered in the middle of the domain:
\begin{equation} \label{2D_lin_adv_IC} u(x,y,0) = \exp\left(-50(x-0.5)^2\right)\exp\left(-50(y - 0.5)^2\right). \end{equation}

As described in section \ref{sec:2D} we will now solve the two-dimensional kinetic equation \eqref{eq:kin_eq_multidim} with the Maxwellian given by \eqref{eq:2D_maxwellian}. The velocity discretization is now determined by the orthogonal velocity method as follows: we fix $R=S=1$ and calculate $v_{\textrm{max}}$ by choosing the (integer) lower bound from expression \eqref{eq:2D_v_max},
\[ v_{\textrm{max}} = \left\lceil \sqrt{\frac{12R^2(a^2+b^2)}{(R+1)(2R+1)}}\;\right\rceil. \]
The inner integrator is a space-time discretization of the kinetic equation \eqref{eq:kin_eq_multidim}, in which we take the ENO scheme of order 1, 2 and 3 in space with $\Delta x = \Delta y = 0.04$ and the forward Euler scheme in time with $\delta t = \epsi = 10^{-8}$. The outer integrator is the fourth order projective Runge-Kutta (PRK4) method, using $K=2$ and $\Delta t = 0.3\Delta x$.

We compare the obtained numerical results for increasing order in space from $1$ to $3$. This is shown in figure \ref{fig:2D_linadv_PRK4_FE_ENO}.

\begin{figure}
	\begin{center}
		\includegraphics[scale=0.50]{2D_linadv_PRK4_FE_ENO123_dxdy004}
	\end{center}
	\caption{\label{fig:2D_linadv_PRK4_FE_ENO} Linear advection in 2D on $t = 1$ for different orders in space of the ENO scheme.}
\end{figure}

\subsection{Dam-break problem in 2D}\label{subsec:shawa_2d}
In dam-break problems, one is interested in the evolution of two (or more) regions of water that are separated by a dam, which then is suddenly removed. Such problems are typically modeled by shallow water equations and can be understood as the equivalent of shock tube applications in gas dynamics (see section \ref{subsec:sod}) \cite{leveque2002finite}. In two space dimensions, the shallow water equations take the following form:
\begin{equation}\label{2D_shawa_eq}
	\begin{cases}
		\partial_t h + \partial_x (h\bar v_x) + \partial_y (h\bar v_y) = 0 \\
		\partial_t (h\bar v_x) + \partial_x (h\bar v_x^2 + \pr) + \partial_y (h\bar v_x \bar v_y) = 0 \\
		\partial_t (h\bar v_y) + \partial_x (h\bar v_x \bar v_y) + \partial_y (h\bar v_y^2 + \pr) = 0, \\
	\end{cases}
\end{equation}
in which $h$ is the depth of the water, ${\boldsymbol{\bar v} = \left(\bar v_x,\bar v_y\right)}$ is the (macroscopic) velocity vector, and $\pr$ is the pressure. All these unknown functions depend on $x, y$ and $t$. The system in \eqref{2D_shawa_eq} contains three partial differential equations for four unknown functions: $h, \bar v_x, \bar v_y$ and $\pr$. Therefore we close the system by the following (hydrostatic) equation of state:
\begin{equation} \label{2D_shawa_eq_of_state} \pr = \frac{1}{2}gh^2. \end{equation}
The shallow water system is considered over the spatial domain $(x,y)\in[-2.5,2.5]^2$ and $t\in[0,1.5]$.

The initial solution consists of a cylindrical basin of water surrounded by a dam (see \cite{leveque2002finite}), given by:
\begin{equation} \label{2D_shawa_eq_IC}
	\begin{aligned}
		h(x,y,0) &= \left\{ \begin{array} {l@{\qquad \quad}l} 2 & x^2+y^2 \le 0.5 \\ 1 & \text{otherwise} \end{array} \right.& \\
		\bar v(x,y,0) &= 0.&
	\end{aligned}
\end{equation}
Furthermore, we impose outflow boundary conditions.

For the relaxation method, we use the kinetic equation \eqref{eq:kin_eq_multidim} together with the Maxwellian given by \eqref{eq:2D_maxwellian}. The velocity discretization is determined by the orthogonal velocity method as follows: we fix $R=S=1$ and calculate $v_{\textrm{max}}$ as follows,
\begin{equation} \label{eq:orthog_vel_vmax} v_{\textrm{max}} = \left\lceil \sqrt{\frac{24R^2}{(R+1)(2R+1)}}\;\right\rceil. \end{equation}
In the projective integration framework, the inner integrator is a space-time discretization of the kinetic equation \eqref{eq:kin_eq_multidim}, in which we take the upwind scheme of order 3 in space with $\Delta x = \Delta y = 0.01$ and the forward Euler scheme in time with $\delta t = \epsi = 10^{-8}$. The outer integrator is the fourth order projective Runge-Kutta (PRK4) method, using $K=2$ and $\Delta t = 0.3\Delta x$. The contour plots of the water depth $h$ can be seen in figure \ref{fig:2D_shawa_PRK4_u3}. As in \cite{leveque2002finite}, we find that at $t=1.5$ the depth of the water near the origin stabilizes around $h\approx0.96$.

\begin{figure}
	\begin{center}
		\includegraphics[scale=0.50]{2D_shawa_PRK4_FE_u3_dxdy001_L250_t150_contours_edited}
	\end{center}
	\caption{\label{fig:2D_shawa_PRK4_u3} Contour plots depicting the evolution of the water depth $h$ on $t=1$ (left) and $t=1.5$ (right) in the cylindrical dam-break problem obtained with PRK4 and FE as inner integrator and the third order upwind scheme. }
\end{figure}

\subsection{Euler equations in 2D}\label{subsec:euler_2d}
When extending the Euler equations from 1D to 2D we end up with a system of four equations since there is conservation of momentum along both the $x$- and $y$-direction. This leads to the following system of equations:
\begin{equation}\label{2D_Euler_eq}
	\begin{cases}
		\partial_t \rho + \partial_x (\rho\bar v_x) + \partial_y (\rho\bar v_y) = 0 \\
		\partial_t (\rho\bar v_x) + \partial_x (\rho\bar v_x^2) + \partial_y (\rho\bar v_x \bar v_y + \pr) = 0 \\
		\partial_t (\rho\bar v_y) + \partial_x (\rho\bar v_y^2) + \partial_y (\rho\bar v_x \bar v_y + \pr) = 0 \\
		\partial_t \E + \partial_x (\E + \pr\bar v_x) + \partial_y (\E + \pr\bar v_y) = 0.
	\end{cases}
\end{equation}
In system \eqref{2D_Euler_eq} the unknown functions $\rho$, ${\boldsymbol{\bar v} = \left(\bar v_x,\bar v_y\right)}$, $\pr$ and $\E$ all depend on $x, y$ and $t$. Similarly to the one-dimensional situation, system \eqref{2D_Euler_eq} now consists of four partial differential equations for five unknown functions: $\rho, \bar v_x, \bar v_y, \pr$ and $\E$. Therefore we close the system by the following equation of state:
\begin{equation} \label{2D_eq_of_state} \pr = (\gamma - 1)\left(\E - \frac{1}{2}\rho||\boldsymbol{\bar v}||^2\right), \end{equation}
with $\gamma = 7/5$. We consider the Euler system over the spatial domain $(x,y)\in[-0.5,0.5]^2$ and $t\in[0,0.18]$.

The extension of the initial discontinuous configuration \eqref{Euler_eq_IC} of Sod's shock test in 1D to 2D is given by (see \cite{aregba2000discrete}):
\begin{equation} \label{2D_Euler_eq_IC}
	\begin{aligned}
		\rho(x,y,0) &= \left\{ \begin{array} {l@{\qquad \quad}l} 0.1 & xy \le 0 \\ 1 & \text{otherwise} \end{array} \right.& \\
		\pr(x,y,0) &= \left\{ \begin{array} {l@{\qquad \quad}l} 0.1 & xy \le 0 \\ 1 & \text{otherwise} \end{array} \right.& \\
		\bar v(x,y,0) &= 0,&
	\end{aligned}
\end{equation}
which is also called a double Sod tube. Furthermore, we impose outflow boundary conditions.

For the relaxation method and velocity discretization, we use the same setting as in section \ref{subsec:shawa_2d}.
In the projective integration framework, the inner integrator is a space-time discretization of the kinetic equation \eqref{eq:kin_eq_multidim}, in which we take the ENO scheme of order 3 in space with $\Delta x = \Delta y = 0.01$ and the forward Euler scheme in time with $\delta t = \epsi = 10^{-8}$. The outer integrator is the fourth order projective Runge-Kutta (PRK4) method, using $K=2$ and $\Delta t = 0.3\Delta x$. The results can be seen in figure \ref{fig:2D_sod_PRK4_ENO3} and correspond to those in \cite{aregba2000discrete}.

\begin{figure}
	\begin{center}
		\includegraphics[scale=0.50]{2D_sod_PRK4_FE_ENO3_dxdy001_t018_contours}
	\end{center}
	\caption{\label{fig:2D_sod_PRK4_ENO3} Evolution of several fluid properties in Sod's shock test in 2D obtained with PRK4 with FE as inner integrator and the third order ENO scheme on $t=0.18$.  }
\end{figure}

\section{Conclusions} \label{sec:conclusions}
We presented a general, high-order, fully explicit, relaxation scheme for systems of nonlinear hyperbolic conservation laws in multiple dimensions, by approximating the nonlinear hyperbolic conservation law by a kinetic equation with BGK source term, which is, in turn, discretized and integrated using a projective integration method. After taking a few small (inner) steps with the direct forward Euler method, an estimate of the time derivative is used in an (outer) Runge-Kutta method of arbitrary order.

Unlike other methods based on relaxation \cite{jin1995relaxation,aregba2000discrete}, the projective integration method does not rely on a splitting technique, but only on an appropriate selection of time steps using a naive explicit discretization method.  Its main advantage is its generality and ease of use: implementing the method for a different system of hyperbolic conservation laws only requires changing the definition of the Maxwellian.

We showed that, with an appropriate choice of inner step size, the time step restriction on the outer time step is similar to the CFL condition for the hyperbolic conservation law. Moreover, the number of inner time steps is also independent of the scaling parameter. We analyzed stability and consistency, and illustrated with numerical results on a set of test problems of varying complexity.

\bibliographystyle{plain} \bibliography{refs}

\end{document}

config: arxiv


TITLE: Kernel of matrix polynomials
QUESTION [0 upvotes]: I am told that $f(A)=g(A)h(A), f(A)=0$ is a polynomial of an $n \times n$ matrix $A$.
I am further told that $g(x)$ and $h(x)$ are coprime. I need to show that $\mathbb{R}^n$ is the direct sum of $\text{ker }g(A)$ and $\text{ker }h(A)$.
My attempt at a solution:
Since $g(x)$ and $h(x)$ are coprime, this means that they have no roots or factors in common, so I can deduce that $\text{ker }g(A) \cap \text{ker }h(A) ={0}$.
I now want to show that if a vector $v \in \mathbb{R}^n$ is not in $\text{ker }g(A)$, then it is in $\text{ker }h(A)$ and vice versa.
Since $f(A)=0$, we have $f(A)v=0 \Rightarrow g(A)h(A)v=0$
Because of commutativity of polynomial multiplication, this also means that $h(A)g(A)v=0$.
If $v \notin \text{ker }g(A)$ then $g(A)v\ne 0$ so $h(A)v=0$. Is this true? Am I on the right track?
Any help would be greatly appreciated!

REPLY [3 votes]: .Use the fact that if $g,h$ are relatively prime polynomials then there exist polynomials $p(x),q(x)$ such that $p(x)g(x) + q(x)h(x) = 1$ for all $x$,(Bezout identity), so in particular, $p(A)g(A) + q(A)h(A) = I$.
Now, with this in mind, let $v \in \Bbb R^n$. Then, we may write $v = p(A)g(A)v + q(A)h(A)v$.
Now, it is clear(by commutativity of the polynomial operators in $A$) that :
$$
h(A)[p(A)g(A)v] = p(A)[g(A)h(A)v] = 0 \\
g(A)[q(A)h(A)v] = q(A)[g(A)h(A)v] = 0
$$
Therefore $v$ is a sum of two elements lying in the kernels of $h(A)$ and $g(A)$ respectively. Adding the disjointness of these subspaces gives the required conclusion.


TITLE: Definition of torsion sheaf on reducible spaces
QUESTION [9 upvotes]: I need to discuss torsion-free sheaves on reduced, but possibly reducible spaces. Here "torsion" means "element is annihilated by a non-zero-divisor". The standard references (EGA, Hartshorne, ...) restrict themselves to normal varieties or  irreducible spaces and do not seem to cover the definition in this level of generality.
I saw some papers and survey articles that use torsion-free sheaves on reducibly spaces, but mostly without any discussion of the definition. Some papers even seem to give the "wrong" definition where "torsion" = "annihilated arbitrary non-zero element", which is probably not what the authors had in mind.
Is anyone aware of a reliable reference for the definition and for basic properties of torsion-free sheaves that I could possibly use?
Thank you in advance!

REPLY [5 votes]: Hi Stefan! Torsion-free (coherent) sheaves are pure sheaves of codimension 0; I suggest that you try "The Geometry of Moduli Spaces of Sheaves", by Daniel Huybrechts and Manfred Lehn, where pure sheaves are discussed in detail.