\section{Error Behaviour using Gauß-Lobatto Nodes}\label{sec:Error_Lobatto}

In this section, Gauß-Lobatto-Legendre nodes will be used in the discretisation,
resulting in diagonal norm SBP operators including the boundary nodes. In this
case, multiplication and restriction to the boundary commute and the interpolated
speed $a(x)$ is automatically continuous.
\revsec{Before starting our investigation in this section,
we will briefly summarize our final results for both cases (Gau\ss-Lobatto-Legendre and Gau\ss-Legendre).}

\begin{Result}
\revsec{$\eta(t)$ is a factor which depends on $\epsilon_1$, the values of $a$
and $a'$. If there exits a positive constant $\delta$, such that the mean value of
 $\eta(t)$ can be bounded from below, then there exists a constant $C$ such
that the errors $\epsilon_1^k(t)$ of \eqref{eq:error_basis} satisfy the inequality
\begin{equation*}
||\epsilon_1(t)||_N\leq \frac{1-\exp(-\delta t)}{\delta} C,
\end{equation*}
in the discrete norm $||\cdot||_N$. The total error is bounded in time.}
\end{Result}

\revsec{
In the following, we will derive the exact conditions when the above inequality
is fulfilled and specify in detail what factors play a key role in the definition of $\eta$ and $\delta$.}
\revPP{
We  outline the steps of our analysis.
All steps of the investigation in sections \ref{sec:Error_Lobatto}
and \ref{sec:Error_Gauss} are almost analogous except
that in step 5 we have to consider the different flux functions
\eqref{eq:Edge_based_cen_flux} - \eqref{eq:unsplit_based_up_flux}
in our investigation\footnote{We have an additional error term in section \ref{sec:Error_Gauss}, but this does not change
the major steps of the study.}.
The main stepts are the following:
\begin{itemize}
 \item[1.] We derive an error equation for $\epsilon_1^k$ of \eqref{eq:error_basis} by inserting the error $E^k$ into
 the continuous equation \eqref{eq:Continuous4} for every element.
 \item[2.] By adding zero in a suitable way, we are able to split the equations into
 a continuous and a discrete part.
 \item[3.] We add both parts for every element and get the error behaviour for the total domain.
 \item[4.] We estimate the continuous terms and get an inequality for the error $\epsilon_1$ in the discrete norms.
 \item[5.] We split the terms with the numerical fluxes. In the Gauss-Legendre case (section \ref{sec:Error_Gauss}),
 we have to be careful with respect to
 the used implementation of the numerical fluxes.
 \item[6.] We estimate the long time error behaviour under some assumptions.
\end{itemize}
}
In the following, the error equation for $\epsilon_1^k=\Ip^N(u^k)-U^k$ will be
derived. Starting by considering Gauß-Lobatto-Legendre nodes in our
semidiscretisation and putting
 $u=\Ip^N(u^k)+\epsilon_p^k$ into \eqref{eq:Continuous4} yields
\begin{align*}
  &\frac{\Delta x_k}{2} \est {\partial_t \Ip^N(u^k), \phi^k} +
   \frac{1}{2}\Bigg( a^k \Ip^N(u^k) \phi^k|_{-1}^1- \est{a^k \Ip^N(u^k), \partial_\xi \phi^k} \\
   &+
   \est{ a^k \partial_\xi \Ip^N(u^k),\phi^k}
   +\est{ \Ip^N(u^k) \partial_\xi a^k,\phi^k}  \Bigg)   \\
   =&- \frac{\Delta x_k}{2} \est {\partial_t \epsilon_p^k, \phi^k} +
   \frac{1}{2} \est{a^k \epsilon_p^k, \partial_\xi \phi^k} -\frac{1}{2}
   \est{ a^k \partial_\xi \epsilon_p^k,\phi^k}
   -\frac{1}{2}\est{ \epsilon_p^k \partial_\xi a^k,\phi^k},
\end{align*}
where $\phi^k \in \P^N$ is a polynomial test function.
For Gauß-Lobatto-Legendre nodes, $a^k\epsilon_p^k=0$ at the endpoints,  since the interpolant
is equal to the solution there. Thus, $a^k\epsilon_p^k\phi^k\Big|_{-1}^1=0$.
Using integration-by-parts for  $\est{ a^k\epsilon_p^k, \partial_\xi \phi^k }$ yields
\begin{equation}\label{eq:continuous_error}
 \begin{aligned}
   &\frac{\Delta x_k}{2} \est {\partial_t \Ip^N(u^k), \phi^k} +
   \frac{1}{2}\Bigg( a^k \Ip^N(u^k) \phi^k|_{-1}^1- \est{a^k \Ip^N(u^k), \partial_\xi \phi^k} \\
   &+
   \est{ a^k \partial_\xi \Ip^N(u^k),\phi^k}
   +\est{ \Ip^N(u^k) \partial_\xi a^k,\phi^k}  \Bigg)   \\
   =&- \frac{\Delta x_k}{2} \est {\partial_t \epsilon_p^k, \phi^k} -
   \frac{1}{2} \est{\partial_\xi(a^k \epsilon_p^k),  \phi^k} -\frac{1}{2}
   \est{ a^k \partial_\xi \epsilon_p^k,\phi^k}
   -\frac{1}{2}\est{ \epsilon_p^k \partial_\xi a^k,\phi^k}.
 \end{aligned}
\end{equation}
\revPP{We have to transfer the continuous scalar product from \eqref{eq:continuous_error}
to the discrete ones. Therefore, we are following the ideas from \cite{oeffner2018error},
add zero to the above equation and rearrange these terms.
We will explain this for the first term on the left side of \eqref{eq:continuous_error}
in detail.
The third to fifth terms on the left side are handled analogously and details can be found in the
appendix \ref{sec:Appendix}.
Applying the interpolation operator together with discrete norms results in


\begin{align}
\est{\partial_t \Ip^N(u^k),\phi^k  }=&
\l( \partial_t \vec{\Ip^N(u^k)}, \vec{\phi}^k \r)_N \nonumber \\ &
+\l\{ \est{ \partial_t \Ip^N(u^k), \phi^k }-
  \l( \partial_t \vec{\Ip^N(u^k)}, \vec{\phi}^k \r)_N  \r\}, \label{eq:interpolation}


\end{align}
Now, we are introducing in the factor $Q$ in the above equation which measures
the projection error of a polynomial of degree $N$ to a polynomial of degree $N-1$.
We can rewrite \eqref{eq:interpolation} as
\begin{align*}
 \est{\partial_t \Ip^N(u^k),\phi^k  }=
\l(\vec{\Ip^N(u^k)}, \vec{\phi}^k \r)_N
+\l\{ \est{ Q(u^k), \phi^k }-
  \l(  \vec{Q(u^k)}, \vec{\phi}^k \r)_N  \r\}
\end{align*}
where $Q(u^k):= \partial _t \l( \Ip^N(u^k)-P^m_{N-1} \l( \Ip^N(u^k)\r)\r)$ and
$P^m_{N-1}$ is the orthogonal projection\footnote{More
details can be found in the appendix \eqref{sec:Appendix}.}
of $u$ onto $\P^{N-1}$ using the inner product of $H^m(e^k)$.
We get similar factors ($Q_1-Q_3$) for the other three terms.
Since $u$ and $a$ are bounded, also all of theses values have to be bounded.
}


Finally, the values of the interpolation polynomial at the boundaries of the element ($-1$ and $1$)
can be approximated by a limitation process from
the left side $\Ip^N(u^k)^{-}$ and right side $\Ip^N(u^k)^{+}$.
To simplify the notation, let
\begin{multline}
 \vecfnumk \l(\Ip^N(u^k)^{-},\Ip^N(u^k)^{+}\r)
 \\:=
 \l(\fnum \l(\Ip_R^N(u^{k-1}),\Ip_L^N(u^{k})\r), \fnum \l(\Ip_R^N(u^{k}),\Ip_L^N(u^{k+1})\r)      \r)^T.
\end{multline}
\revPP{ For boundary points included, the interpolation is continuous
(because the exact solution $u$ is continuous)
and all numerical fluxes are exactly the products of the interpolation and the
coefficient values. One obtains
\begin{equation}\label{eq:in_flux}
 \begin{aligned}
   \frac{1}{2} a^k\Ip^N(u^k)\phi^k\bigg|_{-1}^1=& \vec{\phi}^{k,T} \mat{R}[^T] \mat{B}
 \l(\vecfnumk \l(\Ip^N(u^k)^{-},\Ip^N(u^k)^{+}\r)-\frac{1}{2 }\mat{R}\mat{a^k}\vec{u} \r)
 
 
 
 
 .
 \end{aligned}
\end{equation}}


\revPP{Using the above investigation} and putting \eqref{eq:interpolation}--\eqref{eq:in_flux} in
\eqref{eq:continuous_error} results in
\begin{equation}\label{eq:Delta}
\begin{aligned}
 &\quad
 \frac{\Delta x_k}{2} \l( \partial_t \vec{\Ip^N(u^k)}, \vec{\phi}^k \r)_N
 + \vec{\phi}^{k,T} \mat{R}[^T] \mat{B}
 \l(\vecfnumk \l(\Ip^N(u^k)^{-},\Ip^N(u^k)^{+}\r)-\frac{1}{2 }\mat{R}\mat{a^k}\vec{u} \r)
 \\&\quad
 -\frac{1}{2} \l(  \mat{a^k} \vec{\Ip^N(u^k)} , \partial_\xi \vec{\phi}^k\r)_N
 +\frac{1}{2} \l(  \mat{a^k}  \partial_\xi \vec{\Ip^N(u^k)}, \vec{\phi}^k   \r)_N
 +\frac{1}{2} \l(  \mat{\Ip^N(u^k)} \partial_\xi \vec{a}^k,  \vec{\phi}^k   \r)_N
 \\&=
 +\frac{\Delta x_k}{2} \est{T^k(u), \phi^k}
 +\frac{\Delta x_k}{4} \Bigg( (\vec{Q(u^k)}, \vec{\phi}^k)_N
 -(\vec{Q_1(u^k)}, \partial_x \vec{\phi}^k)_N  \\&\quad
 +(\mat{a^k}\vec{Q_2(u^k)}, \vec{\phi}^k)_N
 +(\mat{Q_3(u^k)}\partial_x \vec{a}^k,  \vec{\phi^k}  )_N   \Bigg)
 +\frac{\Delta x_k}{4}  \est{Q_1(u^k), \partial_x \phi^k},
\end{aligned}
\end{equation}
with
\begin{equation}\label{Def_T}
\begin{aligned}
  T^k(u):= &- \Bigg\{ \partial_t \epsilon_p^k+
 \frac{1}{2} \l( \partial_x(a^k \epsilon_p^k)+\epsilon_p^k \partial_x a^k +a^k \partial_x \epsilon_p^k \r)\\
 &+ \frac{1}{2}
 \l( Q(u^k)+ a^k Q_2(u^k) +(Q_3(u^k) \partial_x a^k) \r) \Bigg\}.
\end{aligned}
\end{equation}
Here, in definition \eqref{Def_T} we have again the derivatives in $x$ since we make the term independent
from the transformation. Therefore, we have in \eqref{eq:Delta} a $\frac{\Delta x_k}{2}$
in the $T^k$ terms.

By \eqref{eq:Gauss-esti_element}, the interpolation error $\epsilon_p^k$ converges in $N$ to zero,
if $m>1$ and the Sobolev norm of the solution is uniformly bounded in time\footnote{Therefore, we need the initial and boundary conditions in the
model problem \eqref{eq:Model_problem}. }.
Equation \eqref{eq:semi_disc_ver} is subtracted form \eqref{eq:Delta}
and with $\epsilon_1^k =\Ip^N(u^k)-U^k$ one
obtains
\begin{align*}\label{eq:Delta2}


& \frac{\Delta x_k}{2} \l( \partial_t \vec{\epsilon}_1^k, \vec{\phi}^k \r)_N
  + \vec{\phi}^{k,T} \mat{R}[^T] \mat{B}
  \l(\vecfnumk \l( (\epsilon_1^k )^{-},(\epsilon_1^k )^{+}\r)-\frac{1}{2 }\mat{R}\mat{a^k}\vec{\epsilon}_1^k
  \r) \\
  -&\frac{1}{2} \l(  \mat{a^k} \vec{\epsilon}_1^k , \partial_\xi \vec{\phi}^k\r)_N
  +\frac{1}{2} \l(  \mat{a^k}  \partial_\xi \vec{\epsilon}_1^k, \vec{\phi}^k   \r)_N
  +\frac{1}{2} \l(  \mat{\epsilon}_1^k \partial_\xi \vec{a}^k,  \vec{\phi}^k   \r)_N \\
  =&+\frac{\Delta x_k}{2} \est{T^k(u), \phi^k}
   +\frac{\Delta x_k}{4}  \est{Q_1(u^k), \partial_x \phi^k}
  +\frac{\Delta x_k}{4} \Bigg( (\vec{Q(u^k)}, \vec{\phi}^k)_N \\
  &-(\vec{Q_1(u^k)}, \partial_x \vec{\phi}^k)_N
  +(\mat{a^k}\vec{Q_2(u^k)}, \vec{\phi}^k)_N
  +\l(\mat{Q_3(u^k)}\partial_x \vec{a}^k,  \vec{\phi^k}   \right)_N
\Bigg) .
\end{align*}
Putting $\phi^k=\epsilon_1^k$ results in the energy equation
\begin{equation}\label{eq:energy1}
\begin{aligned}
 &\quad \frac{\Delta x_k}{4} \frac{\d}{\d t} ||\epsilon_1^k||_N^2+
 \vec{\epsilon}_1^{k,T} \mat{R}[^T] \mat{B}
 \l(\vecfnumk \l( (\epsilon_1^k )^{-},(\epsilon_1^k )^{+}\r)-\frac{1}{2 }\mat{R}\mat{a^k}\vec{\epsilon}_1^k \r)
  \\&\quad-\frac{1}{2} \l(  \mat{a^k} \vec{\epsilon}_1^k , \partial_\xi \vec{\epsilon}_1^k \r)_N
 +\frac{1}{2} \l(  \mat{a^k}  \partial_\xi \vec{\epsilon}_1^k, \vec{\epsilon}_1^k  \r)_N
 +\frac{1}{2} \l(  \mat{\epsilon}_1^k \partial_x \vec{a}^k,  \vec{\epsilon}_1^k   \r)_N
 \\&=
 +\frac{\Delta x_k}{2} \est{T^k(u), \epsilon_1^k}
  +\frac{\Delta x_k}{4}  \est{Q_1(u^k), \partial_x \epsilon_1^k}
  \\&\quad
 +\frac{\Delta x_k}{4} \Bigg( (\vec{Q(u^k)}, \vec{\epsilon}_1^k)_N
 -(\vec{Q_1(u^k)}, \partial_x \vec{\epsilon}_1^k)_N
 +(\mat{a^k}\vec{Q_2(u^k)}, \vec{\epsilon}_1^k)_N \\&\quad  +\l(\mat{Q_3(u^k)}\partial_x \vec{a}^k,  \vec{\epsilon}_1^k   \right)_N \Bigg).
\end{aligned}
\end{equation}

\revPP{Since $\mat{M}[^T]=\mat{M}$, we get
\begin{equation}\label{eq:Deleting_a_terms}
\begin{split}
 &\frac{1}{2} \l(\mat{a^k} \vec{\epsilon}_1^k, \partial_\xi \epsilon_1^k \r)_N= \frac{1}{2}\vec{\epsilon}_1^{k,T} \mat{a^k}[^T] \mat{M}\mat{D}\epsilon_1^k,\\
  &\frac{1}{2} \l( \partial_\xi \epsilon_1^k, \mat{a^k} \vec{\epsilon}_1^k\r)_N=
  \frac{1}{2} \vec{\epsilon}_1^{k,T} \mat{D}[^T] \mat{M} \mat{a^k} \epsilon_1^k=\frac{1}{2}\vec{\epsilon}_1^{k,T} \mat{a}^{k,T} \mat{M}\mat{D}\epsilon_1^k,
\end{split}
\end{equation}
and}


one obtains in \eqref{eq:energy1}
\begin{equation}\label{eq:energy2}
\begin{aligned}
 &\quad \frac{\Delta x_k}{4}   \frac{\d}{\d t} ||\epsilon_1^k||_N^2+
 \vec{\epsilon}_1^{k,T} \mat{R}[^T] \mat{B}
 \l(\vecfnumk \l( (\epsilon_1^k )^{-},(\epsilon_1^k )^{+}\r)-\frac{1}{2 }\mat{R}\mat{a^k}\vec{\epsilon}_1^k \r)
 \\& +\frac{\Delta x_k}{4} \l(  \mat{\epsilon}_1^k  \vec{\epsilon}_1^k  ,  \partial_x \vec{a}^k \r)_N=
 \frac{\Delta x_k}{2} \est{T^k(u), \epsilon_1^k}
  +\frac{\Delta x_k}{4}  \est{Q_1(u^k), \partial_x \epsilon_1^k}
 \\&
 +\frac{\Delta x_k}{4} \Bigg( (\vec{Q(u^k)}, \vec{\epsilon}_1^k)_N
 -(\vec{Q_1(u^k)}, \partial_x \vec{\epsilon}_1^k)_N
 +(\mat{a^k}\vec{Q_2(u^k)}, \vec{\epsilon}_1^k)_N \\ &
 +\l(\mat{Q_3(u^k)}\partial_x \vec{a}^k,  \vec{\epsilon}_1^k \right)_N \Bigg).
 \end{aligned}
\end{equation}
Summing this up over all elements
\revPP{ and by defining the numerical flux of the error as
$\vecepsilon:=\vecfnumk \l( \l( \epsilon_1^k\r)^-,\l(\epsilon_1^k\r)^+\r)$, the global
energy of the error is
}


\begin{equation}\label{eq:energy3}
\begin{aligned}
   &\quad \frac{1}{2}   \frac{\d}{\d t} \sum_{k=1}^K \frac{\Delta x_k}{2}||\epsilon_1^k||_N^2+\sum_{k=1}^K
 \vec{\epsilon}_1^{k,T} \mat{R}[^T] \mat{B}
\l( \vecepsilon-\frac{1}{2 }\mat{R}\mat{a^k}\vec{\epsilon}_1^k \r)
\\&\qquad
 +\frac{1}{2} \sum_{k=1}^K \frac{\Delta x_k}{2} \l(  \mat{\epsilon}_1^k  \vec{\epsilon}_1^k  ,  \partial_x \vec{a}^k \r)_N
  \\&=
  \sum_{k=1}^K\frac{\Delta x_k}{2} \est{T^k(u), \epsilon_1^k}
  + \sum_{k=1}^K\frac{\Delta x_k}{4}  \est{Q_1(u^k), \partial_x \epsilon_1^k} - \sum_{k=1}^K\frac{\Delta x_k}{4}
  (\vec{Q_1(u^k)}, \partial_x \vec{\epsilon}_1^k)_N
  \\&\quad
 +  \sum_{k=1}^K\frac{\Delta x_k}{4} \left( (\vec{Q(u^k)}, \vec{\epsilon}_1^k)_N
 +(\mat{a^k}\vec{Q_2(u^k)}, \vec{\epsilon}_1^k)_N
 +\l(\mat{Q_3(u^k)}\partial_x \vec{a}^k,  \vec{\epsilon}_1^k \right)_N      \right).
 \end{aligned}
\end{equation}
The right-hand side of \eqref{eq:energy3} will be estimated using the Cauchy-Schwarz
inequality. \revsec{For example (the others terms are handled similarly),
\begin{align}
  \sum_{k=1}^K\frac{\Delta x_k}{2} \est{T^k(u), \epsilon_1^k} &\leq
  \sqrt{ \sum_{k=1}^K\frac{\Delta x_k}{2} ||T^k(u)||^2 }  \sqrt{ \sum_{k=1}^K\frac{\Delta x_k}{2} ||\epsilon_1^k||^2   } ,
  \label{eq:T_est}\\


   \sum_{k=1}^K\frac{\Delta x_k}{4} \l(\vec{Q_1(u^k)}, \partial_x \vec{\epsilon}_1^k\r)_N
   &\leq  \frac{1}{2} \sqrt{ \sum_{k=1}^K\frac{\Delta x_k}{2} ||\vec{Q_1(u^k)}||_N^2 }
   \sqrt{ \sum_{k=1}^K\frac{\Delta x_k}{2} ||\partial_x \vec{\epsilon}_1^k||_N^2     } , \label{Eq:est_Q_1}


\end{align}}
Using an estimation for the differential operator $\partial_x$ and the fact that $\epsilon_1\in \P^N$, it is
$||\partial_x \vec{\epsilon}_1^k||_N^2 \leq c_1N^2   ||\vec{\epsilon}_1^k||_N^2$ with a positive constant  $c_1$.
This is due to the fact that all norms are equivalent and we can estimate with a
Markov-Bernstein type inequality,
see \cite{govil1999markov}.
 The estimation is used \revsec{for example} in \eqref{Eq:est_Q_1}. 
 An alternative approach would have been
 to use the summation-by-parts property \eqref{eq:SBP} and estimate analogously.

With the global norm over all elements and the equivalence between the continuous and discrete norms, we obtain
\begin{equation}\label{eq:R_absch}
\begin{aligned}
&\frac{1}{2}   \frac{\d}{\d t} \sum_{k=1}^K \frac{\Delta x_k}{2}||\epsilon_1^k||_N^2+\sum_{k=1}^K
 \vec{\epsilon}_1^{k,T} \mat{R}[^T] \mat{B}
\l( \vecepsilon-\frac{1}{2 }\mat{R}\mat{a^k}\vec{\epsilon}_1^k \r)
 \\ & +\frac{1}{2} \sum_{k=1}^K \frac{\Delta x_k}{2} \l(  \mat{\epsilon}_1^k  \vec{\epsilon}_1^k  ,
 \partial_x \vec{a}^k \r)_N
 \leq  \Bigg\{ c||T ||+\frac{cN}{2}||Q_1|| +
 \frac{1}{2}\Bigg( ||Q||_N+ N\tilde{c_1}||Q_1||_N \\
 &+||a Q_2||_N +||Q_3\partial_x a||_N \Bigg)   \Bigg\}
 ||\epsilon_1||_N
 \equiv \hat{\Ep}(t,N) ||\epsilon_1||_N
 \end{aligned}
\end{equation}
Applying the same approach like in \cite{kopriva2017error} and splitting the sum
into three parts (one for the left physical boundary, one for the right physical
boundary and a sum over the internal element endpoints), it is
\begin{align*}
 &\quad \sum\limits_{k=1}^K \vec{ \epsilon}_1^{k,T}
  \mat{R}[^T] \mat{B} \l(\vecepsilon -\frac{1}{2} \mat{R}\mat{a^k}\vec{\epsilon_1}^k \r)
  \\&=
 \sum\limits_{k=1}^K \vec{ \epsilon}_1^{k,T}
  \mat{R}[^T] \mat{B} \Bigg(\vecfnumk \l( \l( \epsilon_1^k\r)^-,\l(\epsilon_1^k\r)^+\r) -
  \frac{1}{2} \mat{R}\mat{a^k} \vec{\epsilon}_1^k \Bigg)
 \\&=
 - \Eps_L^1 \l( f^{\mathrm{num},1}_L -\frac{1}{2} a_L^1\Eps_L^1 \r)+
 \sum\limits_{k=2}^K \l( f^{\mathrm{num},k}_L -\frac{1}{2} a^{k-1}_R\left(\Eps_R^{k-1}+\Eps_L^{k}  \right) \r)
 \\&\qquad \l( \Eps_R^{k-1}-\Eps_L^{k}  \r)
 + \Eps_R^K \l(f^{\mathrm{num},K}_R -\frac{1}{2} a^{K}_R \Eps^K_R \r).
\end{align*}
 Here,
 $\Eps^k_i$ ($i=L,R; \; k=1,\dots, K$) represents the error $\epsilon_1^k$ at the
 the position in the elements, and (to shorten the notation)
 $f^{\mathrm{num},k}_L:=f^{\mathrm{num},k} \l(\Eps^{k-1}_R,\Eps^{k}_L \r)$, $ f^{\mathrm{num},1}_L:=f^{\mathrm{num},1} \l(0,\Eps^{1}_L \r)$
 and $f^{\mathrm{num},K}_R:= f^{\mathrm{num},1} \l(\Eps^K_R,0 \r) $.
 The external states for the physical boundary contributions are zero,
 because $\Ip^N(u^1)=g$ at the left boundary and the external state for $U^1$ is set to $g$.
 At the right boundary, where the upwind numerical flux is used,
 it doesn't matter what the external state is,
 since its coefficient in the numerical flux is zero.
 One gets for the inner element with $\jump{\Eps^k}=\Eps_R^{k-1}-\Eps_L^{k} $
 \begin{equation*}
 \begin{aligned}
  \sum\limits_{k=2}^K \l( f^{\mathrm{num},k}_L -\frac{1}{2} a^{k-1}_R\left(\Eps_R^{k-1}+\Eps_L^{k}  \right)
  \r)
  \l( \Eps_R^{k-1}-\Eps_L^{k}  \r)& = \sum\limits_{k=2}^K \frac{\sigma a^{k-1}_R}{2}
  \l(\jump{\Eps^{k} }\r)^2
  \geq0,\quad \\
 &\text{with}
 \begin{cases}
 \sigma=0 \quad \text{central flux, }\\
 \sigma=1 \quad \text{ upwind flux. }  \end{cases}
 \end{aligned}
 \end{equation*}
\revsec{For the left and right boundaries, it is finally
\begin{align*}
 &\text{left:} &-\Eps_L^1 \l( f^{\mathrm{num},1}_L -
 \frac{1}{2}  a_L^1\Eps_L^1 \r)=
 
 
 \frac{\sigma a_L^1}{2}\l(\Eps_L^1\r)^2,
 \\
 &\text{right:}&\Eps_R^K \l(f^{\mathrm{num},K}_R -
 \frac{1}{2} a_R^K\Eps^K_R \r)=
 
 
 
 \frac{\sigma  a_R^K}{2} \l(\Eps^K_R  \r)^2.
\end{align*}}
Therefore, the energy growth rate is bounded by
\begin{equation}\label{eq:Errorbound}
\begin{aligned}
   &\frac{1}{2} \frac{\d}{\d t} ||\epsilon_1||^2_N +
  \underbrace{\frac{\sigma}{2} \l( a^K_R\l(\Eps^K_R  \r)^2 +a_L^1\l(\Eps_L^1\r)^2  \r) +
  \frac{\sigma}{2} \sum\limits_{k=2}^K a^{k-1}_R\l(\jump{\Eps^{k} }\r)^2 }_{BTs }
  \\ &+
  \underbrace{\frac{1}{2} \sum_{k=1}^K \frac{\Delta x_k}{2} \l(  \mat{\epsilon}_1^k
  \vec{\epsilon}_1^k,  \partial_x \vec{a}^k \r)_N}_{Int_d}
  \leq \Ep(t,N)
  ||\epsilon_1||_N.
\end{aligned}
\end{equation}

It is $BTs\geq0$.
If $Int_d\geq 0$, then \eqref{eq:Errorbound} has the same form as in
 in \cite{ nordstrom2007error}
and one may estimate/bound analogously
to \cite{kopriva2017error, nordstrom2007error} the error
in time.
The $\Ep$ term depends also on $N$,
but this has no influence in the estimation here.
We rewrite \eqref{eq:Errorbound} as
\begin{equation}\label{eq:Errorbound2}
\revsec{\frac{\d}{\d t}} ||\epsilon_1||_N +
\underbrace{ \frac{BTs+ Int_d}{  ||\epsilon_1||_N^2} }_{\eta(t) } ||\epsilon_1||_N  \leq \Ep(t).
\end{equation}
Like it was described in \cite{nordstrom2007error},
it is assumed that the mean value of $\eta(t)$ over any
finite time interval is bounded by a positive constant $\delta_0$
from below. This means that
$\overline{\eta}\geq \delta_0>0$.
Under the assumption for $u$, the right hand side $\Ep(t,N)$ is also bounded
in time and one can put $\max\limits_{s\in [0,\infty)} \Ep(s,N)\leq C_1<\infty$.
Applying these facts in \eqref{eq:Errorbound2} and integrating over time,
the following inequality for the error is obtained
\begin{equation}\label{eq:Errorbound3}
 ||\epsilon_1(t)||_N \leq \frac{1-\exp(-\delta_0 t)}{\delta_0} C_1,
\end{equation}
see \cite[Lemma 2.3]{nordstrom2007error} for details.
\begin{remark}\label{remark_int}
The term $Int_d$ is a crucial factor.
 If $\partial_x \vec{a}^k >0$, one may estimate the left side of \eqref{eq:Errorbound2}
using the minimum of the discrete values of $\vec{a}$.
Then, $Int_d\geq \frac{1}{2} \min\{\partial_x \vec{a}^k\} ||\epsilon_1||_N^2>0$ and
the above assumption
on $\eta$ is inevitably fulfilled.\\
Simultaneously, the term $Int_d$ can also destroy the error boundedness
if the derivatives of $a$ are negative.
It depends then on the sum of $BTs$ and $Int_d$. The upwind fluxes
can therefore rescue the error boundedness \eqref{eq:Errorbound3} whereas applying
the central flux ($\sigma= 0)$ will contribute to an unlimited growth of the error.
We demonstrate this in some examples in section \ref{sec:Num} and \revsec{make a first analytical estimation in \ref{subsec:analytical}.}
\end{remark}
{"config": "arxiv", "file": "1806.02018/4_Error_bound_semidiscrete_Lobatto.tex"}