Boundary treatment of implicit-explicit Runge-Kutta method for hyperbolic
systems with source terms

Weifeng Zhao111Department of Applied Mathematics, University of Science and
Technology Beijing, Beijing 100083, China. E-mail<EMAIL_ADDRESS>Research supported by the National Natural Science Foundation of China (No.
11801030, No. 11861131004). , Juntao Huang222Department of Mathematics,
Michigan State University, East Lansing, MI 48824, USA. E-mail:
<EMAIL_ADDRESS>
Abstract

In this paper, we develop a high order finite difference boundary treatment
method for the implicit-explicit (IMEX) Runge-Kutta (RK) schemes solving
hyperbolic systems with possibly stiff source terms on a Cartesian mesh. The
main challenge is how to obtain the solutions at ghost points resulting from
the wide stencil of the interior high order scheme. We address this problem by
combining the idea of using the RK schemes at the boundary and an inverse Lax-
Wendroff procedure. The former preserves the accuracy of the RK schemes and
the latter guarantees the stability. Our method is different from the widely
used approach for the explicit RK schemes by imposing boundary conditions at
intermediate stages [12], which could not be derived for the IMEX schemes. In
addition, the intermediate boundary conditions are only available for explicit
RK schemes up to third order while our method applies to arbitrary order IMEX
and explicit RK schemes. Moreover, the present boundary treatment method may
be adapted to IMEX RK schemes solving many other partial differential
equations. For a specific third-order IMEX scheme, we demonstrate the good
stability and third-order accuracy of our boundary treatment through both 1D
examples and 2D reactive Euler equations.

Keywords: Boundary treatment; hyperbolic systems with source terms; implicit-
explicit Runge-Kutta method; high order finite difference; inverse Lax-
Wendroff procedure

##  1 Introduction

Implicit-explicit (IMEX) Runge-Kutta (RK) schemes [7] are a powerful tool for
the time discretization of many partial differential equations (PDEs), such as
convection-diffusion equations [28, 21, 29], hyperbolic systems with stiff
sources [22, 23, 11, 8, 9, 10], and the Boltzmann equation [13, 16, 15, 19].
They achieve good stability by implicitly treating the stiff terms in the PDEs
and high order accuracy by using several immediate stages. Though powerful and
efficient, the IMEX RK schemes for the time-dependent PDEs are not without
drawbacks. A long-standing problem is the lack of appropriate boundary
treatment method [7], which is fundamental for practical applications of IMEX
RK schemes.

For the explicit RK schemes solving hyperbolic equations, a popular boundary
treatment method is to impose consistent boundary conditions for each
intermediate stage [12]. However, these boundary conditions are derived only
for explicit RK schemes up to third order. Though it is remedied for fourth
order schemes in [1, 24], the methods therein only apply to one-dimensional
(1D) scalar equations or systems with all characteristics flowing into the
domain at the boundary. Moreover, the intermediate-stage boundary conditions
are not available for IMEX RK schemes since the implicit part would result in
nonlinear equations which cannot be analytically solved in general. On the
other hand, boundary conditions for intermediate stages are derived and
analyzed in [2, 3, 4, 5, 6] for implicit RK schemes solving linear scalar
equations. At this point, we remark that another limitation of all the
boundary treatments mentioned above [12, 1, 24, 2, 3, 4, 5, 6] is that they
only deal with 1D case where the boundary is located at a grid point,
restricting their applications for 2D problems on complex domains.

Perhaps due to the difficulties of extending the above boundary treatments of
explicit and implicit schemes to the IMEX case, there are rare works on that
of IMEX RK schemes. To our knowledge, the only literature on this topic is
[29]. In that reference, a strategy of imposing intermediate boundary
conditions is presented for one-dimensional (1D) linear convection-diffusion
equations. There is a similar restriction that it does not apply to general
nonlinear equations.

In this paper, we propose a boundary treatment method for the IMEX RK schemes
solving hyperbolic systems with possibly stiff source terms. High order finite
difference schemes on a Cartesian mesh are used for the space discretization.
Due to the wide stencil of high order schemes, ghost points are needed near
the boundary where the numerical stencil is partially outside of the
computational domain. Then the main task of boundary treatment is to obtain
the solutions at the ghost points. We address this problem by combining the
idea of using the RK schemes at the boundary and an inverse Lax-Wendroff (ILW)
procedure [25, 26]. The ILW procedure is employed in [25, 26] as a technique
of boundary treatment for hyperbolic systems by using the PDEs to convert
spatial derivatives of inflow boundaries to time derivatives, while those of
the outflow boundaries are approximated with the solutions at interior points.
Nevertheless, the boundary treatment in [25, 26] relies on the intermediate-
stage boundary conditions [12], which are not available for IMEX schemes as
aforementioned. To overcome this difficulty, we propose to use the RK schemes
at the boundary directly, instead of imposing intermediate-stage boundary
conditions. Combining this idea and the ILW procedure, we can obtain the
solutions as well as their derivatives at the boundary and then the solutions
of the ghosts points are computed by the Taylor expansion at the boundary. In
this way, our boundary treatment not only preserves the accuracy of the IMEX
RK schemes but also processes good stability. In addition, it allows the
boundary to be arbitrarily located. While our method applies to general IMEX
RK schemes, we demonstrate its effectiveness with a specific third-order
scheme in [7] by fixing the CFL number as 0.8. The third-order accuracy of the
method is verified for both 1D and 2D cases, and the numerical results of the
2D reactive Euler equations are comparable to those obtained by the third-
order positivity-preserving schemes with reflective boundary conditions [17].

It is interesting to note that for 1D scalar equations we show that our method
applied to third-order IMEX RK scheme can yield the widely used intermediate-
stage boundary conditions in [12], indicating the correctness of our boundary
treatment for IMEX RK schemes. Since the IMEX RK schemes include explicit
schemes as special cases, the present boundary treatment also applies to
explicit RK schemes. It does not rely on intermediate-stage boundary
conditions and thus is expected to be valid for explicit RK schemes of order
higher than three. This will be investigated in our next work. Moreover, our
method may be adapted to IMEX RK schemes solving the other partial
differential equations, e.g. convection-diffusion equations [28, 21, 29] and
the Boltzmann equation [13, 16, 15, 19]. Many unresolved issues for the
boundary treatment of RK schemes can be explored based on the idea of this
work.

This paper is organized as follows. In Section 2, we introduce the IMEX RK
schemes for hyperbolic systems with source terms. In Section 3, we use the 1D
systems to illustrate our idea of boundary treatment. In Section 4, we show
that for 1D scalar equations our boundary treatment can yield the widely used
intermediate-stage boundary conditions. We extend our boundary treatment to
the 2D reactive Euler equations in Section 5. Section 6 is devoted to
numerical validations of our method. Finally, some conclusions and remarks are
given in Section 7.

##  2 IMEX RK schemes for hyperbolic systems with source terms

Consider a one-dimensional nonlinear hyperbolic system with source term

$\partial_{t}U+\partial_{x}F(U)=Q(U),$ (2.1)

where $0\leq x\leq 1$ and $U,F(U),Q(U)\in\mathbb{R}^{M}$. In the domain $0\leq
x\leq 1$, appropriate boundary conditions and initial data $U(0,x)$ are
required to close the above system. We assume that the Jacobian matrix
$F_{U}(U(t,0))$ always has $p$ positive eigenvalues and thus $p$ independent
relations among incoming and outgoing modes are given at the left boundary
$x=0$:

$B(U(t,0),t)=0.$ (2.2)

Similarly, the boundary condition at $x=1$ are also imposed properly.

A general IMEX RK scheme for the hyperbolic system (2.1) is [11]

$\displaystyle U^{(i)}=U^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\partial_{x}F(U^{(j)})+\Delta
t\sum_{j=1}^{i}a_{ij}Q(U^{(j)}),\quad i=1,2,\ldots,s,$ (2.3a) $\displaystyle
U^{n+1}=U^{n}-\Delta
t\sum_{i=1}^{s}\tilde{\omega}_{i}\partial_{x}F(U^{(i)})+\Delta
t\sum_{i=1}^{s}\omega_{i}Q(U^{(i)}).$ (2.3b)

Here the $s\times s$ matrices $\tilde{A}=(\tilde{a}_{ij})$ and $A=(a_{ij})$
satisfy $\tilde{a}_{ij}=0$ for $j\geq i$ and $a_{ij}=0$ for $j>i$ so that the
scheme is explicit for the convection part and implicit for the source term.
Along with the coefficient vectors
$\tilde{w}=(\tilde{w}_{1},\tilde{w}_{2},\ldots,\tilde{w}_{s})^{T}$ and
$w=(w_{1},w_{2},\ldots,w_{s})^{T}$, the scheme can be represented by a double
Butcher tableau

$\tilde{c}$ |  $\tilde{A}$
---|---
|  $\tilde{w}^{T}$
$c$ |  $A$
---|---
|  $w^{T}$

where $\tilde{c}=(\tilde{c}_{1},\tilde{c}_{2},\ldots,\tilde{c}_{s})^{T}$,
$c=(c_{1},c_{2},\ldots,c_{s})^{T}$ are defined as

$\tilde{c}_{i}=\sum_{j=1}^{i-1}\tilde{a}_{ij},\quad
c_{i}=\sum_{j=1}^{i}a_{ij}.$ (2.4)

For the space discretization in (2.3), we use the finite difference WENO
scheme on a Cartesian mesh [18], which needs a stencil of several points in
each direction. As a consequence, ghost points are needed near the boundary
$\partial\Omega$ where the numerical stencil is partially outside of $\Omega$.
For explicit RK schemes, solutions at ghost points can be determined as in
[25, 26] with the aid of an ILW procedure and the boundary condition of each
stage constructed in [12]. However, this approach does not apply to IMEX
schemes since the boundary conditions of intermediate stages could not be
analytically derived for general nonlinear systems. Thus boundary treatment of
IMEX RK schemes remains a challenge and we address this problem in the
following.

##  3 Boundary treatment for 1D system

In this section, we use the one-dimensional system (2.1) to illustrate our
idea of boundary treatment for the IMEX RK scheme (2.3), and then extend it to
the two-dimensional case in Section 5. We employ the third-order IMEX RK
scheme for concreteness, while our method applies to general IMEX schemes.
Additionally, we use the third-order finite difference WENO scheme in space
[20] and take $\Delta t=\mathcal{O}(\Delta x)$.

###  3.1 Computation of solutions at ghost points

We focus on the left boundary $x=0$ of the problem (2.1) and the method can be
similarly applied to the right boundary. Following the notations in [26], we
discretize the interval $(0,1)$ by a uniform mesh

$\frac{\Delta x}{2}=x_{0}<x_{1}<\cdots<x_{N}=1-\frac{\Delta x}{2}$ (3.1)

and set $x_{j},j=-1,-2$ as two ghost points near the left boundary $x=0$ (note
that the boundary indeed can be arbitrarily located). Denote by $U_{j}^{n}$
the numerical solution of $U$ at position $x_{j}$ and time $t_{n}$. Assume
that the interior solutions $U_{j},j=0,1,2,\ldots,N$ have been updated from
time level $t_{n-1}$ to time level $t_{n}$. As in [26], we use a third-order
Taylor approximation to construct the values at the ghost points $x_{j}$ for
$j=-1,-2$, namely,

$U_{j}^{n}=\sum_{k=0}^{2}\frac{x_{j}^{k}}{k!}U^{n,(k)},\quad j=-1,-2.$ (3.2)

Here $U^{n,(k)}$ denotes a $(3-k)$-th order approximation of the spatial
derivative at the boundary point $\frac{\partial^{k}U}{\partial
x^{k}}\big{|}_{x=0,t=t_{n}}$. With this formula, $U_{j}^{n}$ at the ghost
points can be obtained once $U^{n,(k)},k=0,1,2$ are provided, which is well
treated in [25, 26] (see also the next subsection for details). Then we can
obtain $U_{j}^{(1)}$ at the interior points by the RK scheme (2.3a). The next
task is to compute $U_{j}^{(1)},j=-1,-2$. We propose an approach to do this as
follows.

Similar to the above procedure, we compute $U_{j}^{(1)}$ with the third-order
Taylor expansion at the boundary point $x_{b}=0$:

$U_{j}^{(1)}=\sum_{k=0}^{2}\frac{x_{j}^{k}}{k!}U^{(1),(k)},\quad j=-1,-2,$
(3.3)

where $U^{(1),(k)}$ denotes a $(3-k)$-th order approximation of the spatial
derivative $\frac{\partial^{k}U^{(1)}}{\partial x^{k}}\big{|}_{x=0}$. Then we
turn to compute $U^{(1),(k)}$ for $k=0,1,2$. To this end, we apply the first
stage of the RK solver (2.3a) for $U^{(1)}$ at the boundary point $x_{b}=0$:

$U^{(1)}(x_{b})=U^{n}(x_{b})-\tilde{a}_{10}\Delta
t\partial_{x}F(U^{n}(x_{b}))+\Delta ta_{10}Q(U^{(1)}(x_{b})).$ (3.4)

Notice that $U^{n}(x_{b})=U^{n,(0)}+\mathcal{O}(\Delta x^{3})$ and
$\partial_{x}F(U^{n}(x_{b}))=F_{U}(U^{n}(x_{b}))\partial_{x}U^{n}(x_{b})=F_{U}(U^{n,(0)})U^{n,(1)}+\mathcal{O}(\Delta
x^{2})$ are already known. Substituting these into (3.4) and solving the
algebra equation, we can obtain an approximation of $U^{(1)}(x_{b})$ and
denote it by $U^{(1),(0)}$. In addition, the error between $U^{(1),(0)}$ and
$U^{(1)}(x_{b})$ defined by (3.4) is $\mathcal{O}(\Delta x^{3})$.

Furthermore, taking derivatives with respect to $x$ on both sides of (3.4)
yields

$\frac{\partial U^{(1)}}{\partial
x}\big{|}_{x=x_{b}}=\partial_{x}U^{n}(x_{b})-\tilde{a}_{10}\Delta
t\partial_{xx}F(U^{n}(x_{b}))+\Delta
ta_{10}Q_{U}(U^{(1)}(x_{b}))\frac{\partial U^{(1)}}{\partial
x}\big{|}_{x=x_{b}}.$ (3.5)

Here $\partial_{x}U^{n}(x_{b})=U^{n,(1)}+\mathcal{O}(\Delta x^{2})$ and
$\partial_{xx}F(U^{n}(x_{b}))=F_{UU}(U^{n}(x_{b}))\partial_{x}U^{n}(x_{b})\partial_{x}U^{n}(x_{b})+F_{U}(U^{n}(x_{b}))\partial_{xx}U^{n}(x_{b})=F_{UU}(U^{n,(0)})U^{n,(1)}U^{n,(1)}+F_{U}(U^{n,(0)})U^{n,(2)}+\mathcal{O}(\Delta
x)$. With these approximations, $\frac{\partial U^{(1)}}{\partial
x}\big{|}_{x=x_{b}}$ can be solved out from (3.5) and the resulting solution
$U^{(1),(1)}$ is a second-order approximation of $\frac{\partial
U^{(1)}}{\partial x}\big{|}_{x=x_{b}}$.

By taking higher order derivatives on both sides of (3.4), one can also
compute $\frac{\partial^{k}U^{(1)}}{\partial x^{k}}\big{|}_{x=x_{b}}$ for
$k=2$. However, this procedure is quite complicated as it involves Jacobian of
Jacobian. Here we simply approximate $\frac{\partial^{k}U^{(1)}}{\partial
x^{k}}\big{|}_{x=x_{b}}$ for $k=2$ by using the $(3-k)$-th order WENO type
extrapolation with $U^{(1)}$ at interior points in [26], which will be
introduced in subsection 3.3. In this way, we obtain $U^{(1),(k)}$ for
$k=0,1,2$ with accuracy of order $(3-k)$. Then $U_{j}^{(1)}$ for $j=-1,-2$ can
be computed by (3.3) with third-order accuracy. Having $U^{(1)}$ at the ghost
points, we can then evolve from $U^{(1)}$ to $U^{(2)}$ using the interior
difference scheme.

Repeating the same procedure for each $U^{(i)},2\leq i\leq s$, we can compute
the solution at the ghost points in the $i$-th intermediate stage. Then we can
update the solution at all interior points in the $(i+1)$-th stage. Finally,
we obtain $U^{n+1}$, i.e., the solution at the end of a complete RK cycle.

###  3.2 Computation of $U^{n,(k)}$ at the boundary

For the sake of completeness, we provide the method in [26] for computing
$U^{n,(k)},k=0,1,2$, i.e., the $(3-k)$-th order approximation of
$\frac{\partial^{k}U^{n}}{\partial x^{k}}$ at the boundary $x=0$.

####  3.2.1 $k=0$

We first do a local characteristic decomposition to determine the inflow and
outflow boundary conditions as in [26]. Denote the Jacobian matrix of the flux
evaluated at $x=x_{0}$ by

$A(U_{0}^{n})=\partial_{U}F(U)\big{|}_{U=U_{0}^{n}}$

and assume that it has $p$ positive eigenvalues
$\lambda_{1},\lambda_{2},\ldots,\lambda_{p}$ and $(M-p)$ negative eigenvalues
$\lambda_{p+1},\lambda_{p+2},\ldots,\lambda_{M}$ with $l_{1},l_{2},...l_{p}$
and $l_{p+1},l_{p+2},...l_{M}$ the corresponding left eigenvectors,
respectively. Define by $V_{j,m}$ the $m$-th component of local characteristic
variable at grid point $x_{j},j=0,1,2$, i.e.,

$V_{j,m}=l_{m}U_{j}^{n},\quad m=1,2,\ldots,M,\quad j=0,1,2.$ (3.6)

We extrapolate the outgoing characteristic variable
$V_{j,m},m=p+1,p+2,\ldots,M$ to the boundary $x_{b}$ with the WENO type
extrapolation [26] (see also subsection 3.3), and denote the extrapolated
$k$-th order derivative by

$V_{x_{b},m}^{*(k)},\quad k=0,1,2.$ (3.7)

With $V_{x_{b},m}^{*(0)}$ and the boundary condition (2.2), $U^{n,(0)}$ at the
boundary can be solved out from the following equation

$\begin{split}&l_{m}U^{n,(0)}=V_{x_{b},m}^{*(0)},\quad m=p+1,p+2,\ldots,M,\\\
&B(U^{n,(0)},t_{n})=0.\end{split}$ (3.8)

####  3.2.2 $k=1$

Having $U^{n,(0)}$, we proceed to compute $U^{n,(1)}$ with the ILW procedure
proposed in [25, 26]. To do this, we take derivative with respect to $t$ for
the boundary condition (2.3)

$B_{U}(U^{n,(0)},t_{n})\partial_{t}U^{n,(0)}+B_{t}(U^{n,(0)},t_{n})=0,$

which can be written as

$B_{U}(U^{n,(0)},t_{n})\partial_{t}U^{n,(0)}=g(U^{n,(0)},t_{n})$

with $g(U^{n,(0)},t_{n}):=-B_{t}(U^{n,(0)},t_{n})$. On the other hand,
multiplying the equation (2.1) with $B_{U}(U^{n,(0)},t_{n})$ from the left
yields

$B_{U}(U^{n,(0)},t_{n})\partial_{t}U^{n,(0)}+B_{U}(U^{n,(0)},t_{n})A(U^{n,(0)})U^{n,(1)}=B_{U}(U^{n,(0)},t_{n})Q(U^{n,(0)}).$

With the above two equations and $V_{x_{b},m}^{*(1)}$ obtained by the WENO
type extrapolation, $U^{n,(1)}$ can be solved out from

$\displaystyle l_{m}U^{n,(1)}=V_{x_{b},m}^{*(1)},\quad m=p+1,p+2,\ldots,M,$
(3.9a) $\displaystyle
B_{U}(U^{n,(0)},t_{n})A(U^{n,(0)})U^{n,(1)}=B_{U}(U^{n,(0)},t_{n})Q(U^{n,(0)})-g(U^{n,(0)},t_{n}).$
(3.9b)

This is the ILW procedure in [25, 26].

####  3.2.3 $k=2$

Following [26], we approximate $\frac{\partial^{2}U^{n}}{\partial
x^{2}}|_{x=0}$ with the WENO type extrapolation. Specifically, since we have
obtained the characteristic variables at grid points near the boundary in
(3.6), the WENO type extrapolation is employed based on these characteristic
variables to compute the second-order derivative $V_{x_{b},m}^{*(2)}$ for each
$m=1,2,\cdots,M$. Then the approximation of $\frac{\partial^{2}U^{n}}{\partial
x^{2}}|_{x=0}$ is given by

$U^{n,(2)}=L^{-1}V_{x_{b},m}^{*(2)},$ (3.10)

where $L$ is local left eigenvector matrix composed of
$l_{1},l_{2},\cdots,l_{M}$. Similarly, $\frac{\partial^{k}U^{n}}{\partial
x^{k}}|_{x=0},k>2$ are also computed with the WENO type extrapolation if
necessary for higher order boundary treatments.

###  3.3 WENO type extrapolation

Finally, we show how to obtain a $(3-k)$-th order approximation of
$\frac{\partial^{k}V}{\partial x^{k}}|_{x=0}$, denoted by $V^{*(k)},k=0,1,2$,
with the stencil $S_{3}=\\{x_{0},x_{1},x_{2}\\}$, which is introduced in [25].
When $V$ is smooth on $S_{3}$, we compute $V^{*(k)}$ as

$V^{*(k)}=\sum_{r=0}^{2}d_{r}\frac{\mbox{d}^{k}p_{r}(x)}{\mbox{d}x^{k}}\big{|}_{x=0},$

where $d_{0}=\Delta x^{2},d_{1}=\Delta x,d_{2}=1-\Delta x-\Delta x^{2}$ and
$p_{r}(x)$ is the $r$-th order Lagrange polynomial constructed with $S_{3}$.
When there is a discontinuity in the stencil $S_{3}$, the WENO type
extrapolation is applied [25]:

$V^{*(k)}=\sum_{r=0}^{2}\omega_{r}\frac{\mbox{d}^{k}p_{r}(x)}{\mbox{d}x^{k}}\big{|}_{x=0},$
(3.11)

where $\omega_{r}$ are nonlinear weights depending on the value of $V_{j}$.
Similar to the classical WENO scheme, the nonlinear weights are given by [25]

$\omega_{r}=\frac{\alpha_{r}}{\sum_{s=0}^{2}\alpha_{s}},\quad\alpha_{r}=\frac{d_{r}}{(\varepsilon+\beta_{r})^{2}}.$

Here $\varepsilon=10^{-6}$ and $\beta_{r}$ are the smoothness indicators
determined by

$\beta_{0}=\Delta
x^{2},\quad\beta_{r}=\sum_{l=1}^{r}\int_{x_{-1}}^{x_{0}}\Delta
x^{2l-1}\left(\frac{\mbox{d}^{l}}{\mbox{d}x^{l}}p_{r}(x)\right)^{2}\mbox{d}x,\quad
r=1,2.$

The explicit expressions of $\beta_{1}$ and $\beta_{2}$ are given in [25].
When $V$ is smooth on $S_{3}$, it is shown in [25] that

$\omega_{0}=\mathcal{O}(\Delta x^{2}),\quad\omega_{1}=\mathcal{O}(\Delta
x),\quad\omega_{2}=1-\omega_{0}-\omega_{1}$

so that the WENO type extrapolation (3.11) is $(3-k)$-th order accurate.

In the end of this section, we summarize our third-order boundary treatment
for the hyperbolic system (2.1) as follows:

Step 1: Compute $U^{n,(k)}$, the $(3-k)$-th order approximation of
$\frac{\partial^{k}U^{n}}{\partial x^{k}}$ at the boundary $x=0$ for $k=0,1,2$
as in subsection 3.2. Then impose $U^{n}_{j}$ at the ghost points $j=-1,-2$
with the Taylor expansion (3.2). This step is the same as that in [26] for
explicit RK schemes. With $U^{n}$ at the ghost points, we can evolve from
$U^{n}$ to $U^{(1)}$ with the interior difference scheme.

Step 2: Compute $U^{(1),(k)}$, the $(3-k)$-th order approximation of
$\frac{\partial^{k}U^{(1)}}{\partial x^{k}}$ at the boundary $x=0$ for
$k=0,1,2$ as in subsection 3.1. Specifically, $U^{(1),(0)}$ and $U^{(1),(1)}$
are solved out from the RK solver (3.4) and (3.5), respectively, and
$U^{(1),(2)}$ is approximated with the WENO type extrapolation in subsection
3.3. Then impose $U^{(1)}_{j}$ at the ghost points $j=-1,-2$ with the Taylor
expansion (3.3). Having $U^{(1)}$ at the ghost points, we can update $U^{(2)}$
at all interior points.

Step 3: Repeat the same procedure as in Step 2 for each $U^{(i)},2\leq i\leq
s$ to compute the solution at the ghost points in the $i$-th intermediate
stage. Then we can update the solution at all interior points in the
$(i+1)$-th stage. Finally, we obtain $U^{n+1}$, i.e., the solution at the end
of a complete RK cycle.

######  Remark 3.1.

Combining the idea of using the RK schemes of intermediate solutions at the
boundary and the ILW procedure is the key of our boundary treatment. The
former preserves the accuracy of the RK method and the latter guarantees the
stability. This is different from the widely used approach for the explicit RK
scheme by imposing boundary conditions for intermediate solutions [12], which
are unavailable for the IMEX schemes. In addition, the intermediate boundary
conditions are only available for explicit RK schemes up to third order while
our method applies to general IMEX RK scheme.

######  Remark 3.2.

Though the idea of using the RK schemes of intermediate solutions has been
stressed in [24], the method therein is only valid for 1D scalar equation or
1D system with all characteristics flowing into the domain at the boundary.
Another limitation of the method in [24] is that the boundary has to be
located at a grid point, which makes it difficult to be extended to 2D case.

######  Remark 3.3.

Since the IMEX RK schemes include explicit schemes as special cases, the
present boundary treatment also applies to explicit RK schemes. It does not
need to impose intermediate boundary conditions and thus is expected to be
valid for explicit RK schemes of order higher than three. This will be
investigated in our next work.

##  4 Analysis of the boundary treatment

In this section, we show that our boundary treatment applied to the third-
order explicit RK scheme solving 1D scalar equations can yield the widely used
boundary conditions for intermediate solutions in [12].

###  4.1 Linear equation

We consider a linear scalar equation

$\displaystyle\partial_{t}u+\partial_{x}u=0,\quad 0\leq x\leq 1,\quad t>0,$
(4.1a) $\displaystyle u(t,0)=g(t)$ (4.1b)

and solve it with the third-order explicit strong-stability-preserving RK
scheme [14]

$\displaystyle u^{(1)}=u^{n}+\Delta t\mathcal{L}(u^{n}),$ (4.2a)
$\displaystyle u^{(2)}=\frac{3}{4}u^{n}+\frac{1}{4}u^{(1)}+\frac{1}{4}\Delta
t\mathcal{L}(u^{(1)}),$ (4.2b) $\displaystyle
u^{n+1}=\frac{1}{3}u^{n}+\frac{2}{3}u^{(2)}+\frac{2}{3}\Delta
t\mathcal{L}(u^{(2)}),$ (4.2c)

where $\mathcal{L}(u)=-\partial_{x}u$. We focus on the boundary treatment at
the left boundary $x_{b}=0$.

According to the boundary treatment in the previous section, we first need to
compute $u^{n}$ and its first- and second-order derivatives at the boundary
point $x=x_{b}$, denoted by $u^{n}_{b}$ and $u^{n,(k)}_{b},k=1,2$,
respectively. For the above scalar equation, $u^{n}_{b}$ is simply given by
the boundary condition:

$u^{n}_{b}=g(t_{n}).$ (4.3)

Then according to subsection 3.2.2, the ILW procedure is employed to compute
the first-order derivative $u^{n,(1)}_{b}$:

$u^{n,(1)}_{b}=\partial_{x}u(t_{n},x_{b})=-\partial_{t}u(t_{n},x_{b})=-g^{\prime}(t_{n}).$
(4.4)

Similarly, the second-order derivative $u^{n,(2)}_{b}$ computed by the ILW
procedure is

$u^{n,(2)}_{b}=g^{\prime\prime}(t_{n}),$ (4.5)

and by the third-order WENO extrapolation is

$u^{n,(2)}_{b}=g^{\prime\prime}(t_{n})+\mathcal{O}(\Delta x).$ (4.6)

Here the solution is assumed to be smooth near the boundary for (4.6).

Next we compute $u^{(1)}_{b}$ and its derivatives $u^{(1),(k)}_{b},k=1,2$. In
our method, the RK scheme for $u^{(1)}$ is applied to determine its value at
the boundary, i.e.,

$u^{(1)}_{b}=u^{n}_{b}-\partial_{x}u^{n}_{b}=g(t_{n})+\Delta
tg^{\prime}(t_{n}),$ (4.7)

where (4.3) and (4.4) have been used. As shown in (3.5), we take derivatives
with respect to $x$ on both sides of (4.2a) and evaluate at $x_{0}$ to obtain

$\partial_{x}u^{(1)}_{b}=\partial_{x}u^{n}_{b}-\Delta
t\partial_{xx}u^{n}_{b},$

and thereby

$u^{(1),(1)}_{b}=-g^{\prime}(t_{n})-\Delta tu^{n,(2)}_{b},$ (4.8)

where (4.4) has been used. Applying the scheme (4.2b) for $u^{(2)}$ at the
boundary $x_{0}$ and using (4.3), (4.7) and (4.8), we arrive at

$\begin{split}u^{(2)}_{b}&=\frac{3}{4}u^{n}_{b}+\frac{1}{4}u^{(1)}_{b}-\frac{1}{4}\Delta
tu^{(1),(1)}_{b},\\\ &=g(t_{n})+\frac{1}{2}\Delta
tg^{\prime}(t_{n})+\frac{1}{4}\Delta t^{2}\Delta tu^{n,(2)}_{b}.\end{split}$
(4.9)

If we use (4.5) for $u^{n,(2)}_{b}$ obtained by the ILW procedure, then

$u^{(2)}_{b}=g(t_{n})+\frac{1}{2}\Delta tg^{\prime}(t_{n})+\frac{1}{4}\Delta
t^{2}g^{\prime\prime}(t_{n}).$ (4.10)

In this case, the relations (4.3), (4.7) and (4.10) are exactly the boundary
conditions imposed at each stage in [12, 25].

If the WENO extrapolation is adopted for $u^{n,(2)}_{b}$, then

$u^{(2)}_{b}=g(t_{n})+\frac{1}{2}\Delta tg^{\prime}(t_{n})+\frac{1}{4}\Delta
t^{2}g^{\prime\prime}(t_{n})+\mathcal{O}(\Delta t^{3}),$ (4.11)

which only differs from that in (4.10) with a third-order term
$\mathcal{O}(\Delta t^{3})$ and thus maintains the third-order accuracy of the
boundary treatment.

###  4.2 Nonlinear equation

We further show that the above analysis also applies to the nonlinear case. To
this end, we use the third-order RK scheme (4.2) for a nonlinear scalar
equation

$\begin{split}&\partial_{t}u+\partial_{x}f(u)=0,\quad 0\leq x\leq 1,\quad
t>0,\\\ &u(t,0)=g(t)\end{split}$

with $f^{\prime}(u(0,t))>0$, and apply our method at the left boundary
$x_{b}=0$.

First, it is trivial that the boundary condition for $u^{n}$ is
$u^{n}_{b}=g(t_{n})$. Using the ILW procedure, we can obtain [25]

$\begin{split}&u^{n,(1)}_{b}=\partial_{x}u(t_{n},x_{b})=-\frac{g^{\prime}(t_{n})}{f^{\prime}(g(t_{n}))}\\\
&u^{n,(2)}_{b}=\partial_{xx}u(t_{n},x_{b})=\frac{f^{\prime}(g(t_{n}))g^{\prime\prime}(t_{n})-2f^{\prime\prime}(g(t_{n}))g^{\prime}(t_{n})^{2}}{f^{\prime}(g(t_{n}))^{3}}.\end{split}$
(4.12)

According to our boundary treatment, the RK scheme (4.2a) for $u^{(1)}$ is
used at the boundary to determine $u^{(1)}_{b}$:

$u^{(1)}_{b}=u^{n}_{b}-\Delta
tf^{\prime}(u^{n}_{b})u_{x}(t_{n},x_{b})=g(t_{n})+\Delta tg^{\prime}(t_{n}).$
(4.13)

Moreover, taking derivatives with respect to $x$ on both sides of the scheme
(4.2a) and evaluating at the boundary point $x=x_{b}$, we have

$\begin{split}u^{(1),(1)}_{b}&=u^{n,(1)}_{b}-\Delta
t\left(f^{\prime\prime}(u^{n}_{b})\partial_{x}u^{2}(t_{n},x_{b})+f^{\prime}(u^{n}_{b})\partial_{xx}u(t_{n},x_{b})\right)\\\
&=-\frac{g^{\prime}(t_{n})}{f^{\prime}(g(t_{n}))}-\Delta
tf^{\prime\prime}(g(t_{n}))\left(\frac{g^{\prime}(t_{n})}{f^{\prime}(g(t_{n}))}\right)^{2}\\\
&\hskip 14.22636pt-\Delta
tf^{\prime}(g(t_{n}))\frac{f^{\prime}(g(t_{n}))g^{\prime\prime}(t_{n})-2f^{\prime\prime}(g(t_{n}))g^{\prime}(t_{n})^{2}}{f^{\prime}(g(t_{n}))^{3}}\\\
&=-\frac{g^{\prime}(t_{n})}{f^{\prime}(g(t_{n}))}-\Delta
t\frac{f^{\prime}(g(t_{n}))g^{\prime\prime}(t_{n})-f^{\prime\prime}(g(t_{n}))g^{\prime}(t_{n})^{2}}{f^{\prime}(g(t_{n}))^{2}}.\end{split}$
(4.14)

Thus $u^{(2)}_{b}$ determined by the scheme (4.2b) is

$\begin{split}u^{(2)}_{b}&=\frac{3}{4}u^{n}_{b}+\frac{1}{4}u^{(1)}_{b}-\frac{1}{4}\Delta
tf^{\prime}(u^{(1)}_{b})u^{(1),(1)}_{b}\\\ &=g(t_{n})+\frac{1}{4}\Delta
tg^{\prime}(t_{n})-\frac{1}{4}\Delta tf^{\prime}(g(t_{n})+\Delta
tg^{\prime}(t_{n}))\\\ &\hskip
14.22636pt\times\left(-\frac{g^{\prime}(t_{n})}{f^{\prime}(g(t_{n}))}-\Delta
t\frac{f^{\prime}(g(t_{n}))g^{\prime\prime}(t_{n})-f^{\prime\prime}(g(t_{n}))g^{\prime}(t_{n})^{2}}{f^{\prime}(g(t_{n}))^{2}}\right)\
\\\ &=g(t_{n})+\frac{1}{4}\Delta t\left(1+\frac{f^{\prime}(g(t_{n})+\Delta
tg^{\prime}(t_{n}))}{f^{\prime}(g(t_{n}))}\right)g^{\prime}(t_{n})\\\ &\hskip
14.22636pt+\frac{1}{4}\Delta t^{2}\frac{f^{\prime}(g(t_{n})+\Delta
tg^{\prime}(t_{n}))\left(f^{\prime}(g(t_{n}))g^{\prime\prime}(t_{n})-f^{\prime\prime}(g(t_{n}))g^{\prime}(t_{n})^{2}\right)}{f^{\prime}(g(t_{n}))^{2}}.\end{split}$
(4.15)

Equations (4.13) and (4.15) are the boundary conditions for $u^{(1)}$ and
$u^{(2)}$, respectively. In particular, (4.13) is the same with (4.7), and
(4.15) degenerates to (4.10) when $f(u)$ is linear. Moreover, for nonlinear
$f(u)$ we conduct the Taylor expansion

$f^{\prime}(g(t_{n})+\Delta tg^{\prime}(t_{n}))=f^{\prime}(g(t_{n}))+\Delta
tg^{\prime}(t_{n})f^{\prime\prime}(g(t_{n}))+\mathcal{O}(\Delta t^{2})$

and then (4.15) reduces to

$\begin{split}u^{(2)}_{b}&=g(t_{n})+\frac{1}{4}\Delta t\left(2+\Delta
t\frac{g^{\prime}(t_{n})f^{\prime\prime}(g(t_{n}))}{f^{\prime}(g(t_{n}))}\right)g^{\prime}(t_{n})\\\
&\hskip 14.22636pt+\frac{1}{4}\Delta
t^{2}\frac{f^{\prime}(g(t_{n}))\left[f^{\prime}(g(t_{n}))g^{\prime\prime}(t_{n})-f^{\prime\prime}(g(t_{n}))g^{\prime}(t_{n})^{2}\right]}{f^{\prime}(g(t_{n}))^{2}}+\mathcal{O}(\Delta
t^{3})\\\ &=g(t_{n})+\frac{1}{2}\Delta tg^{\prime}(t_{n})+\frac{1}{4}\Delta
t^{2}g^{\prime\prime}(t_{n})+\mathcal{O}(\Delta t^{3}).\end{split}$

The above expression of $u^{(2)}_{b}$ only differs from (4.10) with
$\mathcal{O}(\Delta t^{3})$ and thus the overall accuracy is still third
order.

This analysis shows the correctness of our boundary treatment when it is
applied to the explicit third-order RK scheme.

##  5 Boundary treatment for 2D reactive Euler equations

In this section, we extend the above idea of boundary treatment to the 2D
case. For convenience of representation, we concentrate on the 2D reactive
Euler equations

$\partial_{t}U+\partial_{x}F(U)+\partial_{y}G(U)=S(U),\quad(x,y)\in\Omega,$
(5.1)

where

$\begin{split}&U=(\rho,\rho u,\rho v,E,\rho Y)^{T},\\\ &F(U)=(\rho u,\rho
u^{2}+p,\rho uv,(E+p)u,\rho uY)^{T},\\\ &G(U)=(\rho v,\rho uv,\rho
v^{2}+p,(E+p)v,\rho vY)^{T},\\\ &S(U)=(0,0,0,0,\omega)\end{split}$

and $E=\frac{1}{2}\rho(u^{2}+v^{2})+\frac{p}{\gamma-1}+\rho qY$. Here $\rho$
is the density, $u$ and $v$ are the velocities in $x$ and $y$ directions, $E$
is the total energy, $p$ is the pressure, $Y$ is the reactant mass fraction,
$q>0$ is the heat release of reaction and $\gamma$ is the specific heat ratio.
The source term is assumed to be in an Arrhenius form

$\omega=-\tilde{K}\rho Ye^{-\tilde{T}/T},$

where $T=p/\rho$ is the temperature, $\tilde{T}>0$ is the activation constant
temperature and $\tilde{K}>0$ is a constant rate coefficient.

The IMEX RK scheme for the 2D reactive Euler equations (5.1) is

$\displaystyle U^{(i)}=U^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\partial_{x}F(U^{(j)})-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\partial_{y}G(U^{(j)})+\Delta
t\sum_{j=1}^{i}a_{ij}S(U^{(j)}),\quad i=1,2,\ldots,s,$ (5.2a) $\displaystyle
U^{n+1}=U^{n}-\Delta
t\sum_{i=1}^{s}\tilde{\omega}_{i}\partial_{x}F(U^{(i)})-\Delta
t\sum_{i=1}^{s}\tilde{\omega}_{i}\partial_{y}G(U^{(i)})+\Delta
t\sum_{i=1}^{s}\omega_{i}S(U^{(i)}),$ (5.2b)

where the coefficients $\tilde{a}_{ij}$, $a_{ij}$, $\tilde{\omega}_{i}$ and
$\omega_{i}$ are the same with those in (2.3). Assume the solutions of all the
grid points inside $\Omega$ have been updated from time level $t_{n-1}$ to
$t_{n}$. For a ghost point $P$, we find a point
$P_{0}=(x_{0},y_{0})=\bm{x}_{0}$ on the boundary $\partial\Omega$ such that
the normal $\bm{n}(\bm{x}_{0})$ at $P_{0}$ goes through $P$. Here
$\bm{n}(\bm{x}_{0})$ points to the exterior of $\Omega$. We set up a local
coordinate system at $P_{0}$ by

$\left(\begin{array}[]{c}\hat{x}\\\
\hat{y}\end{array}\right)=\left(\begin{array}[]{cc}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{array}\right)\left(\begin{array}[]{c}x\\\
y\end{array}\right)=\mathsf{T}\left(\begin{array}[]{c}x\\\
y\end{array}\right),$ (5.3)

where $\theta$ is the angle between the normal $\bm{n}(\bm{x}_{0})$ and the
$x$-axis, and $\mathsf{T}$ is a rotational matrix. The $\hat{x}$-axis then
points in the same direction as $\bm{n}(\bm{x}_{0})$ and the $\hat{y}$-axis
points in the tangent direction. In this local coordinate system, the reactive
Euler equations (5.1) become

$\hat{U}_{t}+F(\hat{U})_{\hat{x}}+G(\hat{U})_{\hat{y}}=S(\hat{U})$ (5.4)

with

$\hat{U}=(\rho,\rho\hat{u},\rho\hat{v},E,\rho
Y)^{T},\quad\left(\begin{array}[]{c}\hat{u}\\\
\hat{v}\end{array}\right)=\mathsf{T}\left(\begin{array}[]{c}u\\\
v\end{array}\right).$

The IMEX scheme in this coordinate is

$\displaystyle\hat{U}^{(i)}=\hat{U}^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\partial_{\hat{x}}F(\hat{U}^{(j)})-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\partial_{\hat{y}}G(\hat{U}^{(j)})+\Delta
t\sum_{j=1}^{i}a_{ij}S(\hat{U}^{(j)}),\quad i=1,2,\ldots,s,$ (5.5a)
$\displaystyle\hat{U}^{n+1}=\hat{U}^{n}-\Delta
t\sum_{i=1}^{s}\tilde{\omega}_{i}\partial_{\hat{x}}F(\hat{U}^{(i)})-\Delta
t\sum_{i=1}^{s}\tilde{\omega}_{i}\partial_{\hat{y}}G(\hat{U}^{(i)})+\Delta
t\sum_{i=1}^{s}\omega_{i}S(\hat{U}^{(i)}).$ (5.5b)

For a third-order boundary treatment, $\hat{U}^{n}$ at the ghost point $P$ is
imposed by the Taylor expansion

$\hat{U}^{n}(P)=\sum_{k=0}^{2}\frac{\Delta^{k}}{k!}\hat{U}^{n,(k)},$ (5.6)

where $\Delta$ is the $\hat{x}$-coordinate of $P$ and $\hat{U}^{n,(k)}$ is a
$(3-k)$-th order approximation of the normal derivative
$\frac{\partial^{k}\hat{U}}{\partial\hat{x}^{k}}\big{|}_{(x,y)=\bm{x}_{0},t=t_{n}}$.
With $\hat{U}^{n}$ at ghost points, $\hat{U}^{(1)}$ at interior points can be
obtained with the scheme (5.5a). Since the computation of $\hat{U}^{n,(k)}$ is
the same with that in [26] and its basic idea has been illustrated for the 1D
case, we omit the details here.

The next step is to compute $\hat{U}^{(1)}$ at the ghost points. Similar to
the 1D case, we use the Taylor expansion

$\hat{U}^{(1)}(P)=\sum_{k=0}^{2}\frac{\Delta^{k}}{k!}\hat{U}^{(1),(k)}$ (5.7)

with $\hat{U}^{(1),(k)}$ a $(3-k)$-th order approximation of
$\frac{\partial^{k}\hat{U}^{(1)}}{\partial\hat{x}^{k}}\big{|}_{(x,y)=\bm{x}_{0}}$,
and turn to determine $\hat{U}^{(1),(k)}$, $k=0,1,2$. First,
$\hat{U}^{(1),(0)}=\hat{U}^{(1)}(P_{0})$ can be computed with scheme (5.5a),
i.e.,

$\hat{U}^{(1)}(P_{0})=\hat{U}^{n}(P_{0})-\Delta
t\tilde{a}_{11}\partial_{\hat{x}}F(\hat{U}^{n})(P_{0})-\Delta
t\tilde{a}_{11}\partial_{\hat{y}}G(\hat{U}^{n})(P_{0})+\Delta
ta_{11}S(\hat{U}^{(1)})(P_{0}).$ (5.8)

Furthermore, taking derivatives with respect to $\hat{x}$ and $\hat{y}$ on
both sides of the above equation, we obtain

$\begin{split}&\partial_{\hat{x}}\hat{U}^{(1)}(P_{0})=\partial_{\hat{x}}\hat{U}^{n}(P_{0})-\Delta
t\tilde{a}_{11}\partial_{\hat{x}\hat{x}}F(\hat{U}^{n})(P_{0})-\Delta
t\tilde{a}_{11}\partial_{\hat{x}\hat{y}}G(\hat{U}^{n})(P_{0})+\Delta
ta_{11}S_{U}(\hat{U}^{(1)})\partial_{\hat{x}}\hat{U}^{(1)}(P_{0}),\\\
&\partial_{\hat{y}}\hat{U}^{(1)}(P_{0})=\partial_{\hat{y}}\hat{U}^{n}(P_{0})-\Delta
t\tilde{a}_{11}\partial_{\hat{x}\hat{y}}F(\hat{U}^{n})(P_{0})-\Delta
t\tilde{a}_{11}\partial_{\hat{y}\hat{y}}G(\hat{U}^{n})(P_{0})+\Delta
ta_{11}S_{U}(\hat{U}^{(1)})\partial_{\hat{y}}\hat{U}^{(1)}(P_{0}).\end{split}$
(5.9)

Since $\hat{U}^{(1)}(P_{0})$ has already been obtained, and
$\partial_{\hat{x}\hat{x}}F(\hat{U}^{n})(P_{0}),\partial_{\hat{x}\hat{y}}G(\hat{U}^{n})(P_{0}),\partial_{\hat{x}\hat{y}}F(\hat{U}^{n})(P_{0})$
and $\partial_{\hat{y}\hat{y}}G(\hat{U}^{n})(P_{0})$ can be approximated with
$\hat{U}^{n}$ at interior grid points by the 2D WENO type extrapolation [26],
then $\partial_{\hat{x}}\hat{U}^{(1)}(P_{0})$ and
$\partial_{\hat{y}}\hat{U}^{(1)}(P_{0})$ can be solved out from the above
equation. High order derivatives,
$\partial_{\hat{x}\hat{x}}\hat{U}^{(1)}(P_{0})$,
$\partial_{\hat{x}\hat{y}}\hat{U}^{(1)}(P_{0})$ and
$\partial_{\hat{y}\hat{y}}\hat{U}^{(1)}(P_{0})$, are also computed by the 2D
WENO type extrapolation [26]. Then $\hat{U}^{(1)}(P)$ can be computed with
(5.7), and similarly $\hat{U}^{(k)}(P),k\geq 2$ can be obtained. The
velocities at the ghost point can be transformed back to the original
coordinate by

$\left(\begin{array}[]{c}u\\\
v\end{array}\right)=\mathsf{T}^{-1}\left(\begin{array}[]{c}\hat{u}\\\
\hat{v}\end{array}\right)$

while the scalar quantities remain unchanged in the two coordinates. In this
way, both the solutions $U^{n}$ at time $t_{n}$ and intermediate solutions
$U^{(i)}$ at ghost points are determined. Note that the 2D WENO type
extrapolation has been been introduced detailedly in [26] and thus is omitted
here.

##  6 Numerical validations

In this section we conduct several numerical experiments for both 1D and 2D
problems to demonstrate the numerical stability and accuracy of our boundary
treatment. In all the computations, we use the third-order finite difference
WENO scheme [20] with the Lax-Friedrichs flux splitting to form the numerical
fluxes [18] and the third-order IMEX RK scheme proposed in [7] with the double
Butcher tableau

0 |  0 |  0 |  0 |  0 |  0
---|---|---|---|---|---
1/2 |  1/2 |  0 |  0 |  0 |  0
2/3 |  11/18 |  1/18 |  0 |  0 |  0
1/2 |  5/6 |  -5/6 |  1/2 |  0 |  0
1 |  1/4 |  7/4 |  3/4 |  -7/4 |  0
|  1/4 |  7/4 |  3/4 |  -7/4 |  0
0 |  0 |  0 |  0 |  0 |  0
---|---|---|---|---|---
1/2 |  0 |  1/2 |  0 |  0 |  0
2/3 |  0 |  1/6 |  1/2 |  0 |  0
1/2 |  0 |  -1/2 |  1/2 |  1/2 |  0
1 |  0 |  3/2 |  -3/2 |  1/2 |  1/2
|  0 |  3/2 |  -3/2 |  1/2 |  1/2

Additionally, we fix the CFL number to be $0.8$ in all the numerical examples.
For all the 1D examples, the mesh for computation is taken as (3.1), namely,
the boundary is located in the middle of two grid points. Different locations
of the boundary are studied for 2D problems.

###  6.1 1D examples

Example 1. We first consider a 1D scalar hyperbolic equation

$\partial_{t}u+u\partial_{x}u=u^{2}+u,$ (6.1)

in the domain $0\leq x\leq 1$ with analytical solution $u=\exp(t+x)$. Since
$u>0$, the boundary condition only needs to be imposed at $x=0$. Specifically,
the initial and boundary conditions are set as

$u(x,0)=\exp(x),\qquad u(0,t)=\exp(t).$ (6.2)

We compute the $L^{1},L^{2}$ and $L^{\infty}$ errors at the terminal time
$t=1$, and the results are listed in Table 6.1. It can be seen that the
convergence order is around three.

Table 6.1: Example 1: error table for the scalar nonlinear equation (6.1).  $\Delta x$ |  $L^{1}$ error | order |  $L^{2}$ error | order |  $L^{\infty}$ error | order
---|---|---|---|---|---|---
1/20 | 7.62e-4 |  | 1.29e-3 |  | 4.22e-3 |
1/40 | 1.13e-4 | 2.75 | 2.17e-4 | 2.57 | 9.57e-4 | 2.14
1/80 | 1.41e-5 | 3.01 | 2.59e-5 | 3.07 | 1.07e-4 | 3.16
1/160 | 1.50e-6 | 3.23 | 2.42e-6 | 3.42 | 7.20e-6 | 3.89
1/320 | 1.57e-7 | 3.25 | 2.32e-7 | 3.38 | 5.53e-7 | 3.70

Example 2. Next we solve a 1D linear system

$\begin{split}&\partial_{t}u+\partial_{x}v=0\\\
&\partial_{t}v+\partial_{x}u=-\frac{1}{\epsilon}(u+v)\end{split}$ (6.3)

in the domain $0\leq x\leq 1$. The Jacobian matrix has two eigenvalues $\pm 1$
and thus only one boundary condition should be imposed at $x=0$ and also
$x=1$, respectively. We construct two sets of analytical solutions to (6.3):

$u(t,x)=\exp(t+x),\quad v(t,x)=-\exp(t+x)$ (6.4)

and

$u(t,x)=\exp(t+x)+\exp(-\frac{x}{\epsilon}),\qquad v(t,x)=-\exp(t+x).$ (6.5)

There exists no boundary layer for the first solution (6.4) while the boundary
layer occurs at $x=0$ for the second one (6.5) with small $\epsilon>0$. We
test our method for the two cases by taking the initial data to match the
analytical solution. The mesh is given as (3.2) and the solutions of $u$ at
$t=1$ are used to evaluate the errors. In the second case, the three interior
grid points near the left boundary are excluded in the computation of the
errors due to the existence of the boundary layer.

For the solution (6.4), we specify $u$ and $v$ at the left and right
boundaries, respectively, i.e.,

$\begin{split}&u(t,0)=\exp(t)\quad\mbox{at}\quad x=0,\\\
&v(t,1)=-\exp(t+1)\quad\mbox{at}\quad x=1.\end{split}$ (6.6)

Similarly, the boundary condition for (6.5) is given by

$\begin{split}&u(t,0)=\exp(t)+1\quad\mbox{at}\quad x=0,\\\
&v(t,1)=-\exp(t+1)\quad\mbox{at}\quad x=1.\end{split}$ (6.7)

The errors for the first case with solution (6.4) are given in Table 6.2, from
which the desired third-order convergence can be observed. For the second case
with solution (6.5), we take $\epsilon=1$ and $10^{-10}$ to test the
convergence order. The results in Table 6.3 show that the present method of
boundary treatment is third-order accurate for both choices of $\Delta
x\ll\epsilon$ and $\Delta x\gg\epsilon$. Fig. 6.1 shows the solutions near the
left boundary for different $\epsilon$ with $\Delta x=1/80$. It can be seen
that our method behaves well even in the presence of boundary layer.

Table 6.2: Example 2: error table for the linear system (6.3) with solution (6.4).  $\Delta x$ |  $L^{1}$ error | order |  $L^{2}$ error | order |  $L^{\infty}$ error | order
---|---|---|---|---|---|---
1/20 | 9.77e-4 |  | 1.06e-3 |  | 2.37e-3 |
1/40 | 1.08e-4 | 3.03 | 1.25e-4 | 3.08 | 3.20e-4 | 2.89
1/80 | 1.46e-5 | 2.88 | 1.71e-5 | 2.87 | 6.59e-5 | 2.28
1/160 | 1.78e-6 | 3.03 | 1.95e-6 | 3.14 | 5.58e-6 | 3.56
1/320 | 2.35e-7 | 2.93 | 2.52e-7 | 2.95 | 7.24e-7 | 2.95
Table 6.3: Example 2: Error table for the linear system (6.3) with solution (6.5).  |  |  | $\epsilon=1$ |  |  |
---|---|---|---|---|---|---
$\Delta x$ |  $L^{1}$ error | order |  $L^{2}$ error | order |  $L^{\infty}$ error | order
1/20 | 1.01e-3 |  | 1.18e-3 |  | 2.39e-3 |
1/40 | 1.20e-4 | 3.08 | 1.34e-4 | 3.14 | 3.17e-4 | 2.91
1/80 | 1.57e-5 | 2.93 | 1.78e-5 | 2.91 | 6.49e-5 | 2.29
1/160 | 1.90e-6 | 3.05 | 2.02e-6 | 3.14 | 5.52e-6 | 3.56
1/320 | 2.49e-7 | 2.93 | 2.62e-7 | 2.95 | 7.19e-7 | 2.94
|  |  | $\epsilon=10^{-10}$ |  |  |
$\Delta x$ |  $L^{1}$ error | order |  $L^{2}$ error | order |  $L^{\infty}$ error | order
1/20 | 9.66e-4 |  | 1.16e-3 |  | 2.38e-3 |
1/40 | 1.16e-4 | 3.06 | 1.32e-4 | 3.13 | 3.22e-4 | 2.89
1/80 | 1.53e-5 | 2.93 | 1.77e-5 | 2.90 | 6.58e-5 | 2.29
1/160 | 1.84e-6 | 3.04 | 2.01e-6 | 3.14 | 5.72e-6 | 3.52
1/320 | 2.51e-7 | 2.88 | 2.68e-7 | 2.91 | 7.07e-7 | 3.02

Figure 6.1: Example 2: Comparison of the computational results (markers) and
analytical solution (line) of the linear system (6.3) with different
$\epsilon$ and fixed $\Delta x=1/80$.

Example 3. We further construct the following nonlinear system

$\begin{split}&\partial_{t}u+\partial_{x}v=0\\\
&\partial_{t}v+\partial_{x}[u+\frac{1}{2}(u+v)^{2}]=-\frac{1}{\epsilon}[u+v+(u+v)^{2}]\end{split}$
(6.8)

in the domain $0\leq x\leq 1$ to test the accuracy of our method. The
expressions (6.4) and (6.5) are also analytical solutions to (6.8) if
cooperated with appropriate initial and boundary conditions. The Jacobian
matrix of the flux is

$A=\left(\begin{array}[]{cc}0&1\\\ 1+u+v&u+v\end{array}\right)$

and it can be diagonalized as

$LAR=\left(\begin{array}[]{cc}1+u+v&0\\\ 0&-1\end{array}\right),$

where $1+u+v>0$ and

$L=\left(\begin{array}[]{cc}1&1\\\ 1&-\frac{1}{1+u+v}\end{array}\right),\quad
R=L^{-1}=\left(\begin{array}[]{cc}\frac{1}{2+u+v}&\frac{1+u+v}{2+u+v}\\\
\frac{1+u+v}{2+u+v}&-\frac{1+u+v}{2+u+v}\end{array}\right).$

The boundary conditions in (6.6) and (6.7) are employed here.

The errors of the solution (6.4) with boundary condition (6.6) are given in
Table 6.4 and the desired third-order accuracy is achieved again. Similar to
the linear system, we take $\epsilon=1$ and $10^{-10}$ for the solution (6.5)
with boundary condition (6.7) and the errors are shown in Table 6.5. It is
clear that the convergence orders are about 3 for both choices of $\epsilon$.
These results demonstrate the third-order accuracy of our method for nonlinear
systems.

Table 6.4: Example 3: error table for the nonlinear system (6.8) with solution (6.4).  $\Delta x$ |  $L^{1}$ error | order |  $L^{2}$ error | order |  $L^{\infty}$ error | order
---|---|---|---|---|---|---
1/20 | 9.25e-4 |  | 1.09e-3 |  | 2.39e-3 |
1/40 | 1.13e-4 | 3.03 | 1.29e-4 | 3.08 | 3.22e-4 | 2.89
1/80 | 1.52e-5 | 2.89 | 1.76e-5 | 2.87 | 6.61e-5 | 2.28
1/160 | 1.86e-6 | 3.03 | 2.01e-6 | 3.13 | 5.61e-6 | 3.56
1/320 | 2.44e-7 | 2.93 | 2.60e-7 | 2.95 | 7.28e-7 | 2.95
Table 6.5: Example 3: error table for the nonlinear system (6.8) with solution (6.5).  |  |  | $\epsilon=1$ |  |  |
---|---|---|---|---|---|---
$\Delta x$ |  $L^{1}$ error | order |  $L^{2}$ error | order |  $L^{\infty}$ error | order
1/20 | 2.52e-4 |  | 2.99e-4 |  | 4.33e-4 |
1/40 | 4.28e-5 | 2.56 | 4.77e-5 | 2.65 | 6.74e-5 | 2.68
1/80 | 6.09e-6 | 2.81 | 6.58e-6 | 2.86 | 8.98e-6 | 2.91
1/160 | 7.86e-7 | 2.95 | 8.33e-7 | 2.98 | 1.10e-6 | 3.03
1/320 | 9.25e-8 | 3.09 | 9.66e-8 | 3.11 | 1.22e-7 | 3.18
|  |  | $\epsilon=10^{-10}$ |  |  |
$\Delta x$ |  $L^{1}$ error | order |  $L^{2}$ error | order |  $L^{\infty}$ error | order
1/20 | 9.66e-4 |  | 1.16e-3 |  | 2.38e-3 |
1/40 | 1.16e-4 | 3.06 | 1.32e-4 | 3.13 | 3.22e-4 | 2.89
1/80 | 1.53e-5 | 2.93 | 1.77e-5 | 2.90 | 6.58e-5 | 2.29
1/160 | 1.84e-6 | 3.04 | 2.01e-6 | 3.14 | 5.72e-6 | 3.52
1/320 | 2.51e-7 | 2.88 | 2.68e-7 | 2.91 | 7.07e-7 | 3.02

###  6.2 2D reactive Euler equations

All the following examples are devoted to validate our boundary treatment for
the 2D reactive Euler equations (5.1). The parameters are fixed as in [27]:

$\gamma=1.2,\quad q=50,\quad\tilde{T}=50,\quad\tilde{K}=2566.4.$

Example 4. We first verify the accuracy of the method through a problem with
analytical solution

$\rho=1+0.3\sin[2\pi(x+y-t)],\quad u=1,\quad v=0,\quad p=1,\quad Y=0$ (6.9)

in a square domain $[0,1]\times[0,1]$, subject to periodic boundary conditions
at the left and right boundaries and $v=0$ at the top and bottom boundaries.
The latter boundary condition is treated with our method. The mesh for
computation is illustrated in Fig. 6.2. As shown in the figure, the distance
between the bottom boundary and its adjacent gird line is $\eta\Delta x$ and
that at the top boundary is $(1-\eta)\Delta x$. The left and right boundaries
are placed in a similar way. We test the accuracy with $\eta=0.5,0.7$ and 0.9,
and the results are listed in Table 6.6. It can be seen that the influence of
$\eta$ on the errors is slight and the convergence order is about 3 for all
the three choices of $\eta$. These demonstrate the third-order accuracy of our
method for 2D problems.

Figure 6.2: Example 4: Illustration of the mesh. The thin lines are the grid lines and the thick lines are the boundaries of the square domain $[0,1]\times[0,1]$. The distance between the bottom boundary and its adjacent gird line is $\eta\Delta x$ and that at the top boundary is $(1-\eta)\Delta x$. In a similar way the left and right boundaries are placed. Two layers of grid points are required outside each boundary for the third-order WENO scheme. Table 6.6: Example 4: error table for 2D reactive euler equations with solution (6.9). |  |  | $\eta=0.5$ |  |  |
---|---|---|---|---|---|---
$\Delta x$ |  $L^{1}$ error | order |  $L^{2}$ error | order |  $L^{\infty}$ error | order
1/20 | 1.37e-2 |  | 1.64e-2 |  | 3.18e-2 |
1/40 | 4.12e-3 | 1.73 | 5.55e-3 | 1.56 | 1.29e-2 | 1.30
1/80 | 7.78e-4 | 2.40 | 1.27e-3 | 2.13 | 3.75e-3 | 1.78
1/160 | 7.35e-5 | 3.40 | 1.35e-4 | 3.23 | 5.25e-4 | 2.84
|  |  | $\eta=0.7$ |  |  |
$\Delta x$ |  $L^{1}$ error | order |  $L^{2}$ error | order |  $L^{\infty}$ error | order
1/20 | 1.41e-2 |  | 1.64e-2 |  | 3.10e-2 |
1/40 | 4.10e-3 | 1.78 | 5.56e-3 | 1.56 | 1.27e-2 | 1.28
1/80 | 7.79e-4 | 2.40 | 1.27e-3 | 2.13 | 3.76e-3 | 1.76
1/160 | 7.34e-5 | 3.41 | 1.35e-4 | 3.23 | 5.32e-4 | 2.82
|  |  | $\eta=0.9$ |  |  |
$\Delta x$ |  $L^{1}$ error | order |  $L^{2}$ error | order |  $L^{\infty}$ error | order
1/20 | 1.44e-2 |  | 1.66e-2 |  | 3.08e-2 |
1/40 | 4.12e-3 | 1.81 | 5.58e-3 | 1.57 | 1.26e-2 | 1.28
1/80 | 7.70e-4 | 2.42 | 1.27e-3 | 2.14 | 3.69e-3 | 1.78
1/160 | 7.34e-5 | 3.39 | 1.35e-4 | 3.23 | 5.27e-4 | 2.81

Example 5. Next we test our method with the example in [27] containing
discontinuity. The computational domain is $[0,2]\times[0,2]$ and the initial
data is

$(\rho,u,v,p,Y)=\left\\{\begin{array}[]{l}(1,0,0,80,0.2),\quad x^{2}+y^{2}\leq
0.36,\\\ (1,0,0,10,0.8),\quad\mbox{otherwise},\end{array}\right.$ (6.10)

where $p$ and $Y$ are different from that in [27] to avoid possibly negative
pressure or density. The boundary condition for the left and bottom boundaries
are $u=0$ and $v=0$, respectively, and the solutions of the grid points
outside the right and top boundaries are equal to those at the grid points
inside the domain.

The mesh for computation is the same as that illustrated in Fig. 6.2 with a
ratio $\eta$. We compare our results with those of the third-order positivity-
preserving scheme in [17] with reflective boundary conditions (left and
bottom) under the same mesh size $\Delta x=1/40$. The contour plots of $p$ and
$\rho$ at $t=0.16$ are given in Fig. 6.3 and the solutions along the
horizontal diagonal line $y=x$ are plotted in Fig. 6.4. Good agreement between
the two methods can be observed. Additionally, we compare in Fig. 6.5 the
solutions of different ratio $\eta$. It can be seen again that the influence
of $\eta$ is small.

Figure 6.3: Example 5: Contour plots of pressure (top) and density (bottom) of
the present method (left) and the third-order positivity-preserving scheme in
[17] (right) with $\Delta x=1/40$ at $t=0.16$.

Figure 6.4: Example 5: Comparison of the pressure (left) and density (right)
of the present method (circle) and the third-order positivity-preserving
scheme in [17] (line) along the diagonal line $y=x$.

Figure 6.5: Example 5: Comparison of the pressure (left) and density (right)
of the present method with different $\eta$ along the diagonal line $y=x$.

Example 6. We test the detonation diffraction in this example with the same
configuration of boundaries in [27] (see Fig. 6.6 below). The spatial domain
is $[0,5]\times[0,5]$ and the initial condition is

$(\rho,u,v,E,Y)=\left\\{\begin{array}[]{l}(11,6.18,0,970,1),\quad x<0.5,\\\
(5,0,0,400,1),\quad\mbox{otherwise},\end{array}\right.$

which is constructed based on that in [27] to avoid negative pressure or
density. The boundary condition is $\bm{n}\cdot\bm{u}=0$ with $\bm{n}$ being
the normal direction of the boundary, except that at $x=0$,
$(\rho,u,v,E,Y)=(11,6.18,0,970,1)$. The terminal time is $t=0.6$ and the mesh
size is $\Delta x=1/20$. The ratio $\eta$ is fixed as 0.5 for all the
boundaries.

The results of the present boundary treatment and the third-order positivity-
preserving scheme in [17] with reflective boundary conditions are compared in
Figs. 6.6 (contour plot) and 6.7 (solutions along $y=2.5$). It is clear that
good agreements between the two methods are achieved again.

Figure 6.6: Example 6: Contour plots of pressure (top) and density (bottom) of
the present method (left) and the third-order positivity-preserving scheme in
[17] (right) with $\Delta x=1/20$ at $t=0.6$.

Figure 6.7: Example 6: Comparison of the pressure (left) and density (right)
of the present method (circle) and the third-order positivity-preserving
scheme in [17] (line) along the horizontal center line $y=2.5$.

Example 7. Finally, we test our method through an example with multiple
obstacles in [17], which is designed following [27]. In this example, the
spatial domain is $[0,10]\times[0,10]$ and there are two obstacles with
positions $[1,3]\times[0,3]$ and $[5,10]\times[0,5]$, respectively. The
initial condition is set as (6.10) and the boundary condition is
$\bm{n}\cdot\bm{u}=0$ for all the boundaries. We take $\Delta x=1/20$ and
results of the present method and the third-order positivity-preserving scheme
in [17] with reflective boundary conditions are shown in Fig. 6.8 at $t=1$ and
Fig. 6.9 at $t=2$. It can be seen that the results of the two methods are very
similar even with the very coarse mesh, which demonstrate the effectiveness of
our boundary treatment for complex domains.

Figure 6.8: Example 7: Contour plots of pressure (top) and density (bottom) of
the present method (left) and the third-order positivity-preserving scheme in
[17] (right) with $\Delta x=1/20$ at $t=1$.

Figure 6.9: Example 7: Contour plots of pressure (top) and density (bottom) of
the present method (left) and the third-order positivity-preserving scheme in
[17] (right) with $\Delta x=1/20$ at $t=2$.

##  7 Conclusions and remarks

In this paper, we propose a high order finite difference boundary treatment
method for the IMEX RK schemes solving hyperbolic systems with possibly stiff
source terms on a Cartesian mesh. By combining the idea of using the RK
schemes at the boundary and the ILW procedure, our method not only preserves
the accuracy of the RK schemes but also processes good stability. Our method
is different from the widely used approach for the explicit RK schemes by
imposing boundary conditions for intermediate solutions [12], which could not
be derived for the IMEX schemes. In addition, the intermediate boundary
conditions are only available for explicit RK schemes up to third order while
our method applies to arbitrary order IMEX and explicit RK schemes. Moreover,
we show that our method applied to the third-order IMEX RK scheme solving 1D
scalar equations can yield the widely used intermediate-stage boundary
conditions in [12]. Finally, both 1D examples and 2D reactive Euler equations
are used to demonstrate the good stability and third-order accuracy of our
boundary treatment for a specific third-order IMEX scheme.

Since the IMEX RK schemes include explicit schemes as special cases, the
present boundary treatment also applies to explicit RK schemes. It does not
rely on intermediate-stage boundary conditions and thus is expected to be
valid for explicit RK schemes of order higher than three. This will be
investigated in our next work. Furthermore, our method may be adapted to IMEX
RK schemes solving the other partial differential equations, e.g. convection-
diffusion equations [28, 21, 29] and the Boltzmann equation [13, 16, 15, 19].
Many unresolved issues for the boundary treatment of RK schemes can be
explored based on the idea of this work.

## References

  * [1] S. Abarbanel, D. Gottlieb, and M. H. Carpenter.  On the removal of boundary errors caused by Runge–Kutta integration of nonlinear partial differential equations.  SIAM Journal on Scientific Computing, 17(3):777–782, 1996.
  * [2] I. Alonso-Mallo.  Runge-Kutta methods without order reduction for linear initial boundary value problems.  Numerische Mathematik, 91(4):577–603, 2002.
  * [3] I. Alonso-Mallo and B. Cano.  Spectral/Rosenbrock discretizations without order reduction for linear parabolic problems.  Applied Numerical Mathematics, 41(2):247–268, 2002.
  * [4] I. Alonso-Mallo, B. Cano, and M. Moreta.  Order reduction and how to avoid it when explicit Runge–Kutta–Nyström methods are used to solve linear partial differential equations.  Journal of Computational and Applied Mathematics, 176(2):293–318, 2005.
  * [5] I. Alonso-Mallo, B. Cano, and N. Reguera.  Avoiding order reduction when integrating linear initial boundary value problems with lawson methods.  IMA Journal of Numerical Analysis, 37(4):2091–2119, 2016.
  * [6] I. Alonso-Mallo, B. Cano, and N. Reguera.  Avoiding order reduction when integrating linear initial boundary value problems with exponential splitting methods.  IMA Journal of Numerical Analysis, 38(3):1294–1323, 2017.
  * [7] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri.  Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations.  Applied Numerical Mathematics, 25(2-3):151–167, 1997.
  * [8] S. Boscarino, L. Pareschi, and G. Russo.  Implicit-explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit.  SIAM Journal on Scientific Computing, 35(1):A22–A51, 2013.
  * [9] S. Boscarino, L. Pareschi, and G. Russo.  A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscale relaxation.  SIAM Journal on Numerical Analysis, 55(4):2085–2109, 2017.
  * [10] S. Boscarino, J.-M. Qiu, and G. Russo.  Implicit-explicit integral deferred correction methods for stiff problems.  SIAM Journal on Scientific Computing, 40(2):A787–A816, 2018.
  * [11] S. Boscarino and G. Russo.  On a class of uniformly accurate IMEX Runge–Kutta schemes and applications to hyperbolic systems with relaxation.  SIAM Journal on Scientific Computing, 31(3):1926–1945, 2009.
  * [12] M. H. Carpenter, D. Gottlieb, S. Abarbanel, and W.-S. Don.  The theoretical accuracy of runge-kutta time discretizations for the initial boundary value problem: A study of the boundary error.  SIAM Journal on Scientific Computing, 16(6):1241–1252, 1995.
  * [13] F. Filbet and C. Yang.  An inverse Lax-Wendroff method for boundary conditions applied to Boltzmann type models.  Journal of Computational Physics, 245(19):43–61, 2013.
  * [14] S. Gottlieb, C.-W. Shu, and E. Tadmor.  Strong stability-preserving high-order time discretization methods.  SIAM review, 43(1):89–112, 2001.
  * [15] J. Hu, R. Shu, and X. Zhang.  Asymptotic-preserving and positivity-preserving implicit-explicit schemes for the stiff BGK equation.  SIAM Journal on Numerical Analysis, 56(2):942–973, 2018.
  * [16] J. Hu and X. Zhang.  On a class of implicit-explicit Runge-Kutta schemes for stiff kinetic equations preserving the Navier-Stokes limit.  Journal of Scientific Computing, 73(2-3):797–818, 2017.
  * [17] J. Huang, W. Zhao, and C.-W. Shu.  A third-order unconditionally positivity-preserving scheme for production–destruction equations with applications to non-equilibrium flows.  Journal of Scientific Computing, 79(2):1015–1056, 2019.
  * [18] G.-S. Jiang and C.-W. Shu.  Efficient implementation of weighted ENO schemes.  Journal of Computational Physics, 126(1):202–228, 1996.
  * [19] S. Jin, H. Lu, and L. Pareschi.  Efficient stochastic asymptotic-preserving implicit-explicit methods for transport equations with diffusive scalings and random inputs.  SIAM Journal on Scientific Computing, 40(2):A671–A696, 2018.
  * [20] X.-D. Liu, S. Osher, and T. Chan.  Weighted Essentially Non-oscillatory Schemes.  Journal of Computational Physics, 115(1):200–212, 1994.
  * [21] J. Lu, J. Fang, S. Tan, C.-W. Shu, and M. Zhang.  Inverse Lax–Wendroff procedure for numerical boundary conditions of convection–diffusion equations.  Journal of Computational Physics, 317:276–300, 2016.
  * [22] L. Pareschi and G. Russo.  Implicit-explicit Runge–Kutta schemes for stiff systems of differential equations.  Recent Trends in Numerical Analysis, 3:269–289, 2000.
  * [23] L. Pareschi and G. Russo.  Implicit–explicit Runge–Kutta schemes and applications to hyperbolic systems with relaxation.  Journal of Scientific Computing, 25(1):129–155, 2005.
  * [24] D. Pathria.  The correct formulation of intermediate boundary conditions for Runge–Kutta time integration of initial boundary value problems.  SIAM Journal on Scientific Computing, 18(5):1255–1266, 1997.
  * [25] S. Tan and C.-W. Shu.  Inverse Lax-Wendroff procedure for numerical boundary conditions of conservation laws.  Journal of Computational Physics, 229(21):8144–8166, 2010.
  * [26] S. Tan, C. Wang, C.-W. Shu, and J. Ning.  Efficient implementation of high order inverse Lax-Wendroff boundary treatment for conservation laws.  Journal of Computational Physics, 231(6):2510–2527, 2012.
  * [27] C. Wang, X. Zhang, C.-W. Shu, and J. Ning.  Robust high order discontinuous Galerkin schemes for two-dimensional gaseous detonations.  Journal of Computational Physics, 231(2):653–665, 2012.
  * [28] H. Wang, C.-W. Shu, and Q. Zhang.  Stability and error estimates of local discontinuous Galerkin methods with implicit-explicit time-marching for advection-diffusion problems.  SIAM Journal on Numerical Analysis, 53(1):206–227, 2015.
  * [29] H. Wang, Q. Zhang, and C.-W. Shu.  Third order implicit-explicit Runge-Kutta local discontinuous Galerkin methods with suitable boundary treatment for convection-diffusion problems with Dirichlet boundary conditions.  Journal of Computational and Applied Mathematics, 342:164–179, 2018\.