\documentclass[11pt,reqno]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage[colorlinks]{hyperref}
\usepackage{mathrsfs,mathdots}

\usepackage{titlesec}
\titlespacing{\section}{0pt}{2.5ex}{1.5ex}
\titlespacing{\subsection}{0pt}{1.5ex}{1ex}
\titlespacing{\subsubsection}{0pt}{1.5ex}{1ex}
\titleformat{\section}{\large\bfseries\centering}{\thesection}{1em}{}
\titleformat{\subsection}[runin]{\bfseries}{\thesubsection.}{0.5em}{}[.\mbox{\ }]
\titleformat{\subsubsection}[runin]{\bfseries}{\thesubsubsection.}{0.4em}{}[.\mbox{\ }]

%%%%%---------- ORCiD
\newcommand\orcidicon[1]{\href{https://orcid.org/#1}{\includegraphics[scale=0.02]{orcid.pdf}}}

\def\udcs{519.63} %Here the author places classificators of the paper according to Russian classification system
\def\mscs{65M20} %Here the author places  classificators according to the AMS classification list.
\setcounter{page}{1}

\newtheorem{lemma}{Lemma}
\newtheorem{theorem}{Theorem}
\newtheorem{definition}{Definition}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}
\newtheorem{remark}{Remark}
\newtheorem{procedure}{Procedure}

\hypersetup{linkcolor=blue,filecolor=black,urlcolor=blue, citecolor=blue}

\def\logo{\textcolor{blue}{{\bf\Huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}}

\def\semrtop
     {
  \vbox{
     \noindent\logo \hspace{-3mm}  \parbox{12cm}{     \begin{center}
     {\LARGE \textsf{СИБИРСКИЕ ЭЛЕКТРОННЫЕ}} \\[2mm]
     {\LARGE \textsf{МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ}} \\[2mm]
     {\Large \textsf{Siberian Electronic Mathematical Reports}} \\[1mm]
     {\Large\tt{\textcolor{blue}{http://semr.math.nsc.ru}}}\\[1mm]
     \raisebox{1ex}{ISSN 1813-3304 }
     \vspace*{-3mm}
%     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
%     {\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
}
     \vspace{-1mm}
     \noindent
      \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa!} \\
     \hline\hline
     \end{tabular}

     \vspace{3mm}
     {\flushleft\it Vol. 20, No. 2, pp. 144--145 (2023) \hspace{\fill}{\rm\small УДК \udcs \hspace{-4mm}}}
     \newline  %Volume, pages in Russian, filled by Editorial board
     	{\rm\small \textcolor{blue}{https://doi.org/10.33048/semi.2023.20.???}}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaa\;}{\rm\small MSC\ \ \mscs }
  }%\vbox
}


\begin{document}
\renewcommand{\refname}{References}
\renewcommand{\proofname}{Proof.}
\renewcommand{\figurename}{Fig.}

\thispagestyle{empty}

\title[SBP-SAT schemes]{\large Summation-by-parts schemes for symmetric hyperbolic systems}

\author[A.N. MALYSHEV]{{\bf Alexander Malyshev}\protect\orcidicon{0009-0004-8611-2026}}

\address{Aleksandr Malyshev  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} University of Bergen,
\newline\hphantom{iii} Department of Mathematics,
\newline\hphantom{iii} Postbox 7803,
\newline\hphantom{iii} 5020, Bergen, Norway}%
\email{\textcolor{blue}{alexander.malyshev@uib.no}}%

\thanks{\sc Malyshev, A.N.,
Summation-by-parts schemes for symmetric hyperbolic systems}
\thanks{\copyright \ 2024 Malyshev A.N.}
%\thanks{\rm The work is supported by RFFI (grant 03-01-11111)}
\thanks{\it Received January, 1, 2023, Published December, 31, 2023}%

\semrtop \vspace{1cm}
\maketitle {\small

\vspace{-12pt}
\centerline{{\it Communicated by} {\sc P.P. Petrov} %({\it filled by Editorial board})
}

\bigskip
\bigskip

\begin{quote}
\noindent{\bf Abstract:}
We apply the method of lines to numerically solve general initial-boundary value problems for symmetric hyperbolic systems of linear differential equations with variable coefficients. Semi-discretization of symmetric hyperbolic systems is performed using classical summation-by-parts difference operators. Strictly dissipative boundary conditions are weakly enforced using the so-called simultaneous approximation terms.
All theoretical construc\-tions are provided with full proofs.
The stability of explicit Runge-Kutta methods for semi-bounded operators is proved using
recent results on strong stability for semi-dissipative operators.

\medskip

\noindent{\bf Keywords:} symmetric hyperbolic system, dissipative boundary conditions, summation-by-parts scheme, simultaneous approxima\-tion terms, strong stability of explicit Runge-Kutta methods.
\end{quote}
}

\bigskip


\section{Introduction}

We begin our study of the Summation-By-Parts (SBP) difference operators and enforcement of boundary conditions with the Simultaneous Approximati\-on Terms (SAT) with an example demonstrating how the SBP-SAT method extends energy estimates obtained for partial differential equations to semi-discrete approximations of the differential equations with little effort.

\subsection{The energy identity for the linear advection equation}
Solution $u(t,x)$ of the initial-boundary value problem
\begin{align}
&\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}=f(t,x),\quad 0<x<1,\; 0<t<T, \label{f1}\\
&u(0,x)=u_0(x),\quad u(t,0)=g(t), \label{f2}
\end{align}
satisfies the equations $(u^2)_t=u_tu+uu_t=-u_xu+fu-uu_x+uf$ and
\[
\frac{d}{dt}\int_0^1u^2dx=-\int_0^1(u_xu+uu_x)dx+2\int_0^1fudx.
\]
\emph{Integration by parts} ensures the identity
\begin{equation} \label{f3}
\int_0^1(u_xu+uu_x)dx=u^2(1)-u^2(0),
\end{equation}
which allows us to derive the \emph{energy identity}
\begin{equation} \label{f4}
\frac{d}{dt}\int_0^1u^2(t,x)dx=u^2(t,0)-u^2(t,1)+2\int_0^1f(t,x)u(t,x)dx.
\end{equation}

We split the solution $u(t,x)$ into the sum $u=u_1+u_2$, where $u_1(t,x)$ satisfies (\ref{f1}), (\ref{f2}) with $g(t)\equiv0$ while $u_2(t,x)$ satisfies the system (\ref{f1}), (\ref{f2}) with $f(t,x)\equiv0$ and $u_0(x)\equiv0$. It follows from (\ref{f4}) that
\begin{align*}
\frac{d}{dt}\int_0^1u_1^2(t,x)dx&\leq2\int_0^1f(t,x)u_1(t,x)dx\\
&\leq2\left(\int_0^1f^2(t,x)dx\right)^{\frac{1}{2}}\left(\int_0^1u_1^2(t,x)dx\right)^{\frac{1}{2}}
\end{align*}
and, for $0\leq t\leq T$,
\begin{equation} \label{f5}
\left(\int_0^1u_1^2(t,x)dx\right)^{\frac{1}{2}}
\leq\left(\int_0^1u_1^2(0,x)dx\right)^{\frac{1}{2}}+
t\max_{0\leq s\leq T}\left(\int_0^1f^2(s,x)dx\right)^{\frac{1}{2}}.
\end{equation}
Similarly, (\ref{f4}) implies the estimates
\begin{equation} \label{f6}
\frac{d}{dt}\int_0^1u_2^2(t,x)dx \leq g^2(t),\quad \int_0^1u_2^2(t,x)dx\leq\int_0^tg^2(t)dt.
\end{equation}
Combining the estimates (\ref{f5}) and (\ref{f6}),
we arrive at the \emph{energy estimate} for the solution of (\ref{f1}), (\ref{f2}),
\begin{equation} \label{f7}
\begin{split}
\left(\int_0^1u^2(t,x)dx\right)^{\frac{1}{2}}&\leq \left(\int_0^1u_1^2(t,x)dx\right)^{\frac{1}{2}}+\left(\int_0^1u_2^2(t,x)dx\right)^{\frac{1}{2}}
\\
&\leq \left(\int_0^1u_0^2(x)dx\right)^{\frac{1}{2}}\\
&\quad +\left(\int_0^tg^2(t)dt\right)^{\frac{1}{2}}+
t\max_{0\leq s\leq T}\left(\int_0^1f^2(s,x)dx\right)^{\frac{1}{2}}.
\end{split}
\end{equation}

\subsection{Semi-discretization of the advection equation}
Our goal is to reproduce all the derivations of Section 1.1 for an appropriate discrete analo\-gue of (\ref{f1}). Discretization of (\ref{f1}) is carried out only along the $x$ axis, and here we do not use discretization along the variable $t$. The resulting system of ordinary differential equations in $t$ is referred to as a \emph{semi-discretization}. The combination of semi-discretization with a subsequent time integration by Runge-Kutta methods or by multi-step methods is often called the \emph{method of lines}.
In what follows, we consider semi-discrete functions $u(t,x_j)=u_j(t)$ defined at the grid nodes $x_j=jh$, $j=0,1,\ldots,N$, $h=1/N$.

Finite difference approximations of the first derivative of a discrete function $V=\begin{bmatrix}v_0&v_1&\ldots&v_N\end{bmatrix}^T$ are conveniently represented
by the product $\mathsf{D}V$ with a square matrix $\mathsf{D}$, which is often called the \emph{differentiation matrix}.
We choose the differentiation matrix in the factorized form $\mathsf{D}=P^{-1}Q$ with a diagonal matrix $P$ and an ``almost'' skew-symmetric matrix $Q$ as
\begin{equation} \label{f8}
P = \!\begin{bmatrix}\frac{h}{2}\\&h\\&&h\\&&&\ddots\\&&&&h\\&&&&&\frac{h}{2}\end{bmatrix}\!\!,\,
Q = \!\begin{bmatrix}-1/2&1/2\\-1/2&0&1/2\\&-1/2&0&1/2\\&&\ddots&\ddots&\ddots\\
&&&-1/2&0&1/2\\&&&&-1/2&1/2\end{bmatrix}\!\!.
\end{equation}
The semi-discretization of (\ref{f1}) with initial value $u_0$ but without boundary conditions is
as follows:
\begin{equation} \label{f9}
\frac{d}{dt}V(t)+P^{-1}QV(t)=f_h(t),\quad V(0)=(u_0)_h.
\end{equation}
The subscript $h$ indicates the discretization of the functions $f$ and $u_0$ of the continuous argument $x$ at the grid points $x_j$.

The boundary condition $u(t,0)=g(t)$ is imposed by means of the so-called \emph{Simultaneous Approximation Terms (SAT)} method \cite{3}, which enforces boundary conditions in the semi-discretized differential equations as penalty terms. The final semi-discretization is
\begin{equation} \label{f10}
\frac{d}{dt}U(t)+P^{-1}QU(t)=f_h(t)+\mathbb{SAT},\quad U(0)=(u_0)_h
\end{equation}
with the SAT term, depending on a suitable scalar $\sigma$,
\begin{equation} \label{f11}
\mathbb{SAT}=-\frac{\sigma}{2}P^{-1}
\begin{bmatrix}U_0(t)-g(t)\\0\\\vdots\\0\end{bmatrix}.
\end{equation}

\subsection{Convergence analysis for the semi-discrete approximation
of the advection equation}

Let us denote the semi-discretization of the exact solution $u(t,x)$ of (\ref{f1}), (\ref{f2}),
and its derivative $u_x(t,x)$ by the vector functions $u_h$ and $(u_x)_h$ respectively.
Then, obviously, $(u_t)_h=(u_h)_t$ and $(u_h)_t+(u_x)_h=f_h$.
The standard argument based on Taylor expansions proves the consistency equation
\begin{equation} \label{f12}
P^{-1}Qu_h=(u_x)_h+T_h
\end{equation}
with the truncation error
\begin{equation} \label{f13}
T_h=[O(h),O(h^2),\ldots,O(h^2),O(h)]^T.
\end{equation}
As a result, we arrive at the equation
\begin{equation} \label{f14}
\frac{d}{dt}u_h+P^{-1}Qu_h=f_h+T_h.
\end{equation}
In addition, $(u_h)_0=g$ owing to the boundary condition.

The error function $\mathscr{E}(t)=U(t)-u_h(t)=\begin{bmatrix}\mathscr{E}_0&\mathscr{E}_1&\ldots&\mathscr{E}_N\end{bmatrix}^T$ for $U(t)$ satisfying the semi-discretization (\ref{f10}), (\ref{f11}) solves the system
\begin{equation} \label{f15}
\frac{d}{dt}\mathscr{E}+P^{-1}Q\mathscr{E}=
-T_h-\frac{\sigma}{2}P^{-1}\begin{bmatrix}\mathscr{E}_0 & 0 &\ldots & 0\end{bmatrix}^T.
\end{equation}
In the formulas below we use the weighted inner product $(U,V)_P=U^*PV$, the Euclidean inner product $(U,V)=U^*V$ and the weighted norm $\|U\|_P=\sqrt{(U,U)_P}$. We derive from (\ref{f15}) that
\begin{align*}
\frac{d}{dt}(\mathscr{E},\mathscr{E})_P&=(\mathscr{E}_t,\mathscr{E})_P+(\mathscr{E},\mathscr{E}_t)_P\\
&=-(P^{-1}Q\mathscr{E},\mathscr{E})_P-(\mathscr{E},P^{-1}Q\mathscr{E})_P-2(T_h,\mathscr{E})_P-\sigma \mathscr{E}_0^2\\
&=-(Q\mathscr{E},\mathscr{E})-(\mathscr{E},Q\mathscr{E})-2(T_h,\mathscr{E})_P-\sigma \mathscr{E}_0^2\\
&=-((Q+Q^T)\mathscr{E},\mathscr{E})-2(T_h,\mathscr{E})_P-\sigma \mathscr{E}_0^2\\
&=\mathscr{E}_0^2-\mathscr{E}_N^2-2(T_h,\mathscr{E})_P-\sigma \mathscr{E}_0^2\\
&=\mathscr{E}_0^2(1-\sigma)-\mathscr{E}_N^2-2(T_h,\mathscr{E})_P.
\end{align*}
Thus we arrive at the energy equality for the error $\mathscr{E}(t)$,
\begin{equation} \label{f16}
\frac{d}{dt}(\mathscr{E},\mathscr{E})_P=\mathscr{E}_0^2(1-\sigma)-\mathscr{E}_N^2-2(T_h,\mathscr{E})_P.
\end{equation}
It follows from (\ref{f16}) that $\frac{d}{dt}(\mathscr{E},\mathscr{E})_P\leq 2|(T_h,\mathscr{E})_P|$, when $\sigma\geq1$.
Furthermore, $\frac{d}{dt}\|\mathscr{E}\|_P^2\leq 2\|T_h\|_P\|\mathscr{E}\|_P$, $\frac{d}{dt}\|\mathscr{E}\|_P\leq\|T_h\|_P$ and
$\|\mathscr{E}(t)\|_P\leq\|\mathscr{E}(0)\|_P+t\|T_h\|_P$. Note that $\|T_h\|_P=O(h^{3/2})$.
Taking into account $\mathscr{E}(0)=0$, we obtain the estimate of the convergence rate for the scheme (\ref{f10}), (\ref{f11}):
\begin{equation} \label{f17}
\|\mathscr{E}(t)\|_P\leq O(h^{3/2})\mbox{ if } \sigma\geq 1.
\end{equation}
\begin{remark} \label{rem1}
A more refined proof shows that the error $\|\mathscr{E}(t)\|_P$ is bounded from above by $O(h^2)$;
see e.g. \cite{26}.
\end{remark}

\subsection{Summation-By-Parts (SBP) difference operator}

The key com\-ponent in deriving the energy identity (\ref{f4}) is the integration-by-parts rule (\ref{f3}). A convenient discrete analogue of (\ref{f3}) would allow us to convert the derivations of the energy identities for partial differential equations into the semi-discrete case with relative ease. For a differentiation matrix $\mathsf{D}$ and discrete functions $V=\begin{bmatrix}v_0&v_1&\ldots&v_N\end{bmatrix}^T$,
the most convenient discrete analogue of (\ref{f3}) might be the identity $(\mathsf{D}V,V)+(V,\mathsf{D}V)=(\mathsf{E}V,V)$, where $\mathsf{E}=\mbox{diag}[-1,0,0,\ldots,0,1]$ and $(\cdot,\cdot)$ are suitable inner products. This approach indeed works, and we describe it in detail in the present paper.

A differentiation matrix for approximation of the first derivative is sought in the form $\mathsf{D}=P^{-1}Q$, where the matrix $P$ is symmetric positive definite and the matrix $Q$
satisfies the SBP property
\begin{equation} \label{f18}
Q+Q^T=\mathsf{E}=\mbox{diag}[-1,0,0,\ldots,0,1].
\end{equation}
Such a differentiation matrix $\mathsf{D}$ is referred to as the \emph{Summation-By-Parts (SBP)} operator.
The corresponding discrete analogue of the integration-by-parts has the form
\begin{equation} \label{f19}
(P^{-1}QU,V)_P+(U,P^{-1}QV)_P=(\mathsf{E}U,V),
\end{equation}
where $(\cdot,\cdot)$ is the Euclidean inner product and $(\cdot,\cdot)_P$ is the energy inner product $(U,V)_P=(PU,V)$.

\begin{remark} \label{rem2}
Suppose that $a=x_0<x_1<\cdots < x_{s_x}=b$ is an arbitrary grid and consider the vectors $\mathbf{x}^j=\begin{bmatrix}x_0^j&x_1^j&\ldots&x_{s_x}^j\end{bmatrix}^T$. The system of equations $P^{-1}Q\mathbf{x}^j=j\mathbf{x}^{j-1}$ for $0\leq j\leq\tau$ defines accuracy of order $\tau$ for
the SBP operator $P^{-1}Q$. It is equivalent to $Q\mathbf{x}^j=jP\mathbf{x}^{j-1}$.
Multiplying the latter equation from the left by $(\mathbf{x}^i)^T$ leads to the system $(\mathbf{x}^i)^TQ\mathbf{x}^j=j(\mathbf{x}^i)^TP\mathbf{x}^{j-1}$ for $0\leq i,j\leq\tau$.
Now we derive that
\[
(\mathbf{x}^i)^TE\mathbf{x}^j=(\mathbf{x}^i)^TQ\mathbf{x}^j+(\mathbf{x}^i)^TQ^T\mathbf{x}^j=
j(\mathbf{x}^i)^TP\mathbf{x}^{j-1}+i(\mathbf{x}^j)^TP\mathbf{x}^{i-1}.
\]
Let $P=\mbox{\rm diag}\,[p_0,p_1,\ldots,p_{s_x}]$. We arrive at the identity
\begin{equation}\label{f20}
\sum_{k=0}^{s_x}p_kx_k^{i+j-1}=\frac{b^{i+j}-a^{i+j}}{i+j}
\end{equation}
valid for all $i,j$ such that $0\leq i+j\leq 2\tau$. The identity (\ref{f20}) proves that the diagonals $p_k$ are the weights of a quadrature rule of accuracy $2\tau-1$ on the interval $[a,b]$
at the given nodes $x_k$. More about SBP on general grids is found in \cite{5}.
\end{remark}

We use only SBP operators with diagonal matrices $P$ because those with full matrices $P$ are not
applicable to differential operators with variable coefficients. Note that the diagonal matrix $P$
in the SBP operator for the first derivative on an equidistant grid with step size $h$
is equal to a constant diagonal matrix multiplied by $h$.

SBP operators on equidistant grids are called the classical SBP operators.
They were introduced and constructed in \cite{14,15,21}. More recent descrip\-tions
of classical SBP operators can be found in \cite{23,10}. Modern reviews of the SBP-SAT
theory and applications with extensive bibliography are \cite{25,6}. SBP operators
for general function spaces are proposed in \cite{8}.

\section{Classical SBP operators}

\begin{theorem} \label{th1}
The skew-symmetric difference operator on the equidistant mesh $x_{\nu}=\nu h$, $\nu\in\mathbb{N}$,
\begin{equation}\label{f21}
v(x_n)=\frac{1}{h}\left[\sum_{\nu=1}^s\alpha_{\nu}u(x_{n+\nu})-\sum_{\nu=1}^s\alpha_{\nu}u(x_{n-\nu})\right]
\end{equation}
approximates the first derivative $u'(x_n)$ with order $2s$ if and only if
\begin{align}\label{f22}
\sum_{\nu=1}^s\alpha_{\nu}\nu^{2k+1}=\begin{cases}\frac{1}{2},&k=0,\\
0,&k=1,2,\ldots,s-1.\end{cases}
\end{align}
The coefficients $\alpha_{\nu}$ determined by (\ref{f22}) equal
\begin{equation}\label{f23}
\alpha_{\nu}=\frac{(-1)^{\nu-1}(s!)^2}{\nu(s+\nu)!(s-\nu)!},\quad \nu=1,2,\ldots,s.
\end{equation}
\end{theorem}
\begin{proof}
Using the Taylor expansion $u(x_{n+\nu})=\sum_{j=0}^{\infty}u^{(j)}(x_{n})\frac{(\nu h)^j}{j!}$
we obtain the equalities
\begin{align*}
v(x_n)&=\frac{1}{h}\left[\sum_{\nu=1}^s\alpha_{\nu}\sum_{j=0}^{\infty}u^{(j)}(x_{n})\frac{(\nu h)^j}{j!}-\sum_{\nu=1}^s\alpha_{\nu}\sum_{j=0}^{\infty}u^{(j)}(x_{n})\frac{(-\nu h)^j}{j!}\right]\\
&=\frac{1}{h}\sum_{j=0}^{\infty}u^{(j)}(x_{n})\frac{h^j}{j!}[1-(-1)^j]
\sum_{\nu=1}^s\alpha_{\nu}\nu^j=u^{(1)}(x_{n})+O(h^{2s}),
\end{align*}
where the last equality holds because $v(x_n)$ approximates $u'(x_n)$ with order $2s$. Hence
the coefficients $\alpha_{\nu}$ satisfy the system (\ref{f22}).

The system (\ref{f22}) is the Vandermonde system
\[
V\begin{bmatrix}\alpha_1\\2\alpha_2\\\vdots\\s\alpha_s\end{bmatrix}=
\begin{bmatrix}1&1&\ldots&1\\1^2&2^2&\ldots&s^2\\ \vdots&\vdots&\ddots&\vdots\\
(1^2)^{s-1}&(2^2)^{s-1}&\ldots&(s^2)^{s-1}\end{bmatrix}
\begin{bmatrix}\alpha_1\\2\alpha_2\\\vdots\\s\alpha_s\end{bmatrix}=
\begin{bmatrix}1/2\\0\\\vdots\\0\end{bmatrix}.
\]
Cramer's rule gives the solution $\nu\alpha_{\nu}=\frac{1}{2}(-1)^{\nu-1}\det V_{\nu}/\det V$,
where the corresponding submatrix
\[
V_{\nu}=\begin{bmatrix}1^2&\ldots&(\nu-1)^{2}&(\nu+1)^2&\ldots&s^2\\ \vdots&\cdots&\vdots&\vdots&\cdots&\vdots\\
1^{2(s-1)}&\ldots&(\nu-1)^{2(s-1)}&(\nu+1)^{2(s-1)}&\ldots&s^{2(s-1)}\end{bmatrix}
\]
is again of the Vandermonde type. The formulas of the Vandermonde determi\-nants \cite{12} are
$\det V=\Pi_{1\leq i<j\leq s}(j^2-i^2)$ and
\begin{align*}
\det V_{\nu}&=\Pi_{k=1,k\ne\nu}^sk^2\Pi_{1\leq i<j\leq s,\;i\ne\nu,j\ne\nu}(j^2-i^2)\\
&=\frac{(s!)^2}{\nu^2}\frac{\det V}{\Pi_{1\leq i<\nu}(\nu^2-i^2)
\Pi_{\nu<j\leq s}(j^2-\nu^2)}=2\frac{(s!)^2\det V}{(s-\nu)!(s+\nu)!}.
\end{align*}
\end{proof}

We look for a finite difference operator approximating the first derivative on a semi-infinite mesh
$x_{\nu}$, $\nu=0,1,\ldots$, as a block-partitioned semi-infinite matrix
\begin{equation} \label{f24}
\frac{1}{h}\begin{bmatrix}H^{-1}B & H^{-1}C \\ -C^T & D\end{bmatrix}
\end{equation}
having blocks
\begin{align*}
C&=\begin{bmatrix}C_0&0&\cdots\\C_s&0&\cdots\end{bmatrix}, \quad
C_s=\begin{bmatrix}\alpha_s& 0 &\cdots & \cdots & 0\\
\alpha_{s-1} & \alpha_s & 0 & \cdots & 0\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\vdots& & \ddots& \ddots & 0 \\
\alpha_1 & \cdots & \cdots & \alpha_{s-1} & \alpha_s \end{bmatrix},\\
D&=\begin{bmatrix}0 & \alpha_1& \cdots & \alpha_s & 0 & \cdots &&& 0 & \cdots\\
-\alpha_1 & 0 & \alpha_{1} & \cdots & \alpha_s & 0 & \cdots\\
\vdots&\ddots&\ddots&\ddots&&\ddots&\ddots\\
-\alpha_s & \cdots & -\alpha_1 & 0 & \alpha_1 & \cdots & \alpha_s & \ddots && \vdots \\
0 & -\alpha_s & \cdots & -\alpha_1 & 0 & \alpha_1 & \cdots & \alpha_s& 0 & \cdots \\
\vdots & \ddots & \ddots && \ddots & \ddots &\ddots & & \ddots & \ddots \end{bmatrix}.
\end{align*}
The matrix $C_0$ is the zero $(r-s)\times s$ matrix. The coefficients $\alpha_{\nu}$ in $C_s$ and
in the semi-infinite matrix $D$ are defined in Theorem~\ref{th1}. The block $D$ provides $2s$-order accurate approximation of the first derivative at the inner grid points $x_{\nu}$, $\nu\ge r$.

We assume that $r\geq s$ and define the structure of the real $r\times r$ matrix $B$ by
the following conditions:
\begin{equation} \label{f25}
B=B_1+B_2,\quad B_1=\mbox{diag}(-1/2,0,\ldots,0),\quad B_2^T=-B_2.
\end{equation}
We want to find a symmetric positive definite $r\times r$ matrix $H$ and
an ``almost'' skew-symmetric matrix $B$ of the form (\ref{f25})
such that for arbitrary smooth functions $u(x)$,
\begin{equation} \label{f26}
H^{-1}B\begin{bmatrix}u(x_0)\\u(x_1)\\\vdots\\u(x_{r-1})\end{bmatrix}+
H^{-1}\begin{bmatrix}C_0\\C_s\end{bmatrix}
\begin{bmatrix}u(x_r)\\u(x_{r+1})\\\vdots\\u(x_{r+s-1})\end{bmatrix}=
h\begin{bmatrix}u_x(x_0)\\u_x(x_1)\\\vdots\\u_x(x_{r-1})\end{bmatrix}
+\mathcal{O}(h^{\tau+1}).
\end{equation}
Condition (\ref{f26}) provides approximation of order $\tau$ for the first derivative
at the grid nodes $x_{\nu}$, $0\leq\nu<r$, near the boundary. To satisfy (\ref{f26}),
it suffices to check it only for monomial functions $u(x)=(x-x_r)^j$.
It follows that condition (\ref{f26}) is equivalent to the algebraic system
\begin{equation} \label{f27}
(-1)^jH^{-1}B\mathbf{x}^j+(-1)^jH^{-1}\begin{bmatrix}C_0\\C_s\end{bmatrix}\mathbf{y}^j=
(-1)^{j-1}j\mathbf{x}^{j-1}, \quad j=0,1,\ldots,\tau,
\end{equation}
where we denote
\[
\mathbf{x}^j=\begin{bmatrix}r^j&(r-1)^j&\ldots&1^j\end{bmatrix}^T, \quad
\mathbf{y}^j=\begin{bmatrix}0^j&(-1)^j&\ldots&(1-s)^j\end{bmatrix}^T
\]
and agree on that $0^0=1$ and $\mathbf{x}^{-1}=0$.
Multiplying (\ref{f27}) from the left by $H(-1)^j$ and using (\ref{f25}) gives the system
\[
B_2\mathbf{x}^j+B_1\mathbf{x}^j+\begin{bmatrix}C_0\\C_s\end{bmatrix}\mathbf{y}^j=
-jH\mathbf{x}^{j-1}, \quad j=0,1,\ldots,\tau.
\]
As a result, we arrive at the following equations for $j=0,1,\ldots,\tau$
\begin{equation} \label{f28}
B_2\mathbf{x}^j=\mathbf{g}_j, \mbox{ where }\;
\mathbf{g}_j=-B_1\mathbf{x}^j-\begin{bmatrix}C_0\\C_s\end{bmatrix}\mathbf{y}^j
-jH\mathbf{x}^{j-1},
\end{equation}
or in matrix form
\begin{equation} \label{f29}
B_2X=G,
\end{equation}
where $X=\begin{bmatrix}\mathbf{x}^0&\mathbf{x}^1&\cdots&\mathbf{v}^{\tau}\end{bmatrix}$,
$G=\begin{bmatrix}\mathbf{g}_0&\mathbf{g}_1&\cdots&\mathbf{g}_{\tau}\end{bmatrix}$.

\subsection{System of linear equations defining the diagonal matrix $H$}
\begin{theorem} \label{th2}
Under the restriction $\tau+1\leq r$,
a necessary and sufficient condition for existence of a skew-symmetric solution $B_2$ of (\ref{f29})
is the skew-symmetry of the matrix $X^TG$.
\end{theorem}
\begin{proof}
The restriction $\tau+1\leq r$ is needed for existence of a QR factorization $X=QR$, where $R=\begin{bmatrix}\hat{R}\\0\end{bmatrix}$ and $\hat{R}$ is a square non-singular matrix.
The matrix $X$ is a Vandermonde matrix, therefore, $\hat{R}$ is non-singular.
Multiplying the system (\ref{f29}) by $Q^T$ from the left and by $\hat{R}^{-1}$ from the right
gives the equation
\begin{equation} \label{f30}
(Q^TB_2Q)\begin{bmatrix}I\\0\end{bmatrix}=Q^TG\hat{R}^{-1}.
\end{equation}
Let us partition the skew-symmetric matrix $\Theta=Q^TB_2Q$ as
\[
\Theta=\begin{bmatrix}\Theta_{11}&\Theta_{12}\\\Theta_{21}&\Theta_{22}\end{bmatrix},
\]
where the block $\Theta_{11}$ is of size $(\tau+1)\times(\tau+1)$. The block of columns $0,1,\ldots,\tau$ in $\Theta$ equals
\begin{equation} \label{f31}
\begin{bmatrix}\Theta_{11}\\\Theta_{12}\end{bmatrix}=Q^TG\hat{R}^{-1}
\end{equation}
due to equation (\ref{f30}). The upper square block $\Theta_{11}=\begin{bmatrix}I&0\end{bmatrix}Q^TG\hat{R}^{-1}$
of the skew-symmetric matrix $\Theta$ must be skew-symmetric.
Hence the product $\hat{R}^T\Theta_{11}\hat{R}=X^TG$ is also skew-symmetric,
and the necessity is proved.

Sufficiency is proved by using the QR factorization of $X$ again. The block matrix $\begin{bmatrix}\Theta_{11}\\\Theta_{12}\end{bmatrix}$ is constructed as in (\ref{f31}), where
the skew-symmetry of $\Theta_{11}$ is guaranteed by the skew-symmetry of $X^TG$.
The block $\Theta_{12}$ is constructed skew-symmetrically as $\Theta_{12}=-\Theta_{21}^T$. The remaining block $\Theta_{22}$ is chosen arbitrarily but it must be skew-symmetric. Therefore, $\Theta_{22}$ is not unique if its order is larger than 1.
\end{proof}

Skew-symmetry of $X^TG$ is equivalent to the system of equations
\begin{equation} \label{f32}
(\mathbf{x}^i)^Tg_j+g_i^T\mathbf{x}^j=0, \quad i,j=0,1,\ldots,\tau.
\end{equation}
Due to the structure of the vectors $g_i$ in (\ref{f28}) and structure of $B_1$, we have
\[
(\mathbf{x}^i)^T\mathbf{g}_j=\frac{1}{2}r^{i+j}-
(\mathbf{x}^i)^T\begin{bmatrix}C_0\\C_s\end{bmatrix}\mathbf{y}^j
-j(\mathbf{x}^i)^TH\mathbf{x}^{j-1}.
\]
Thus, the system (\ref{f32}) is equivalent to the system
\begin{equation} \label{f33}
j(\mathbf{x}^i)^TH\mathbf{x}^{j-1}+i(\mathbf{x}^j)^TH\mathbf{x}^{i-1}=
r^{i+j}-\xi_{ij}-\xi_{ji},
\end{equation}
where $\xi_{ij}$ are the entries of the matrix $\Xi$ of order $\tau+1$ defined by the formula
\begin{equation} \label{f34}
\xi_{ij}=\begin{bmatrix}s^i&(s-1)^i&\ldots&1^i\end{bmatrix}
C_s\begin{bmatrix}0^j&(-1)^j&\ldots&(1-s)^j\end{bmatrix}^T.
\end{equation}

Recall that we are interested only in the diagonal matrices $H$. Let us denote the diagonal entries of $H$
by the column vector $\mathbf{h}$. In this case,
$(\mathbf{x}^i)^TH\mathbf{x}^{j-1}=(\mathbf{x}^{i+j-1})^T\mathbf{h}$ and
$(\mathbf{x}^j)^TH\mathbf{x}^{i-1}=(\mathbf{x}^{i+j-1})^T\mathbf{h}$ in (\ref{f33}).
The following linear system defining the desired vector $\mathbf{h}$ is obtained from (\ref{f33})
and is equivalent to (\ref{f32}):
\begin{equation} \label{f35}
(i+j)(\mathbf{x}^T)^{i+j-1}\mathbf{h}=r^{i+j}-\xi_{ij}-\xi_{ji},\quad i,j=0,1,\ldots,\tau.
\end{equation}

The system (\ref{f35}) is overdetermined. The left-hand side of (\ref{f35}) is a function of $i+j$, so the right-hand side of (\ref{f35}) must be the same. In other words, the sums $\xi_{ij}+\xi_{ji}$ must be equal for the same $i+j$. A matrix is said to have the Hankel structure if the entries on each anti-diagonal are equal. Thus, for (\ref{f35}) to be consistent, it is necessary that the matrix $\Xi+\Xi^T$ be a Hankel matrix. The equation for $i=j=0$ in (\ref{f35}) is $0=1-2\xi_{00}$. Since this equation does not define $\mathbf{h}$, the equality $\xi_{00}=1/2$ is also a consistency condition.

The number of anti-diagonals in $\Xi+\Xi^T$, excluding the one with $i+j=0$, is $2\tau$. It is equal to the number of independent equations in (\ref{f35}) provided that all consistency conditions are satisfied.
The set of independent equations obviously has a Vandermonde structure and determines $r$ entries of the vector $\mathbf{h}$. Therefore, the maximum $\tau$ in (\ref{f26}) satisfies the equality $2\tau=r$, i.e. the maximum approximation order of the constructed SBP operator is $\tau=s$ provided that the elements of the vector $\mathbf{h}$ satisfying (\ref{f35}) are positive.

\subsection{$\Xi+\Xi^T$ is a Hankel matrix}

For convenience, we introduce the following function of $j$ and $\sigma$, $j\leq\sigma$,
\begin{equation} \label{f36}
J_{\sigma,j}=\xi_{\sigma-j,j}+\xi_{j,\sigma-j},
\end{equation}
and prove that $J_{\sigma,j}$ is independent of $j$. Formula (\ref{f34}) allows us to derive
the equality
\begin{equation} \label{f37}
\xi_{ij}=\sum_{\nu=1}^s\alpha_{\nu}\sum_{\mu=1}^{\nu}\mu^i(\mu-\nu)^j.
\end{equation}
A simple rearrangement of (\ref{f37}) gives the representations
\begin{equation} \label{f38}
\xi_{ji}=\sum_{\nu=1}^s\alpha_{\nu}\sum_{\mu=1}^{\nu}\mu^j(\mu-\nu)^i=
(-1)^{i+j}\sum_{\nu=1}^s\alpha_{\nu}\sum_{\mu=0}^{\nu-1}\mu^i(\mu-\nu)^j
\end{equation}
and
\begin{equation} \label{f39}
\xi_{ij}+\xi_{ji}=\sum_{\nu=1}^s\alpha_{\nu}\left[\sum_{\mu=1}^{\nu}\mu^i(\mu-\nu)^j+
(-1)^{\sigma}\sum_{\mu=0}^{\nu-1}\mu^i(\mu-\nu)^j\right].
\end{equation}
Hence
\begin{equation} \label{f40}
J_{\sigma,j}=\sum_{\nu=1}^s\alpha_{\nu}\left[\sum_{\mu=1}^{\nu}\mu^{\sigma-j}(\mu-\nu)^j+
(-1)^{\sigma}\sum_{\mu=0}^{\nu-1}\mu^{\sigma-j}(\mu-\nu)^j\right].
\end{equation}

\subsubsection{$\sigma$ is odd}

Formula (\ref{f40}) reduces to the equality
\[
J_{\sigma,j}=\sum_{\nu=1}^s\alpha_{\nu}[\nu^i0^j-0^i(-\nu)^j],\quad i=\sigma-j.
\]
Therefore,
\[
J_{\sigma,j}=\begin{cases}\hspace{1em} 0,&i>0,\;j>0,\\
-\sum_{\nu=1}^s\alpha_{\nu}(-\nu)^j,&i=0,\;j>0,\\
\hspace{1em}\sum_{\nu=1}^s\alpha_{\nu}\nu^i,&i>0,\;j=0.\end{cases}
\]
Owing to (\ref{f22}),
\begin{equation} \label{f41}
J_{\sigma,j}=\begin{cases}\frac{1}{2},& \sigma=1,\\
0,&\sigma=3,5,\ldots,2s-1.\end{cases}
\end{equation}

\subsubsection{$\sigma$ is even}

First we consider $\sigma=0$. In this case, $i=j=0$ and
\[
\xi_{ij}+\xi_{ji}=2\sum_{\nu=1}^s\alpha_{\nu}\sum_{\mu=0}^{\nu-1}\mu^0(\mu-\nu)^0=
2\sum_{\nu=1}^s\alpha_{\nu}\nu=1.
\]
Thus, we have proved the consistency condition $\xi_{00}=1/2$.

Let us rewrite (\ref{f40}) in the form
\begin{equation} \label{f42}
J_{\sigma,j}=\sum_{\nu=1}^s\alpha_{\nu}N_{\sigma.j}(\nu),
\end{equation}
where
\begin{equation} \label{f43}
N_{\sigma,j}(\nu)=\sum_{\mu=0}^{\nu-1}\mu^{\sigma-j}(\mu-\nu)^j+
\sum_{\mu=1}^{\nu}\mu^{\sigma-j}(\mu-\nu)^j.
\end{equation}
It is easy to verify that the polynomial $N_{\sigma,j}(\nu)$ of degree $j$
satisfies the recursion
\begin{equation} \label{f44}
N_{\sigma,j}(\nu)=N_{\sigma,j-1}(\nu)-\nu N_{\sigma-1,j-1}(\nu).
\end{equation}
The recursion is solved by the binomial formula
\begin{equation} \label{f45}
N_{\sigma,j}(\nu)=\sum_{n=0}^{j}\binom{j}{n}N_{\sigma-n,0}(\nu)(-\nu)^n.
\end{equation}
\begin{theorem} \label{th3}
\[
N_{\sigma,0}(\nu)=2\begin{cases}
\frac{\nu^{\sigma+1}}{\sigma+1}+\frac{1}{2}B_2\binom{\sigma}{1}\nu^{\sigma-1}+
\frac{1}{4}B_4\binom{\sigma}{3}\nu^{\sigma-3}+\cdots+B_{\sigma}\nu,&\sigma\mbox{ even}\\[1ex]
\frac{\nu^{\sigma+1}}{\sigma+1}+\frac{1}{2}B_2\binom{\sigma}{1}\nu^{\sigma-1}+
\frac{1}{4}B_4\binom{\sigma}{3}\nu^{\sigma-3}+\cdots+\frac{\sigma}{2}B_{\sigma-1}\nu^2,&\sigma\mbox{ odd,}
\end{cases}
\]
where the coefficients $B_{\sigma}$ are Bernoulli's numbers.
\end{theorem}
\begin{proof}
We set $j=0$ in (\ref{f43}) and obtain that
\[
N_{\sigma,0}(\nu)=\sum_{\mu=0}^{\nu-1}\mu^{\sigma}+
\sum_{\mu=1}^{\nu}\mu^{\sigma}=2\sum_{\mu=1}^{\nu}\mu^{\sigma}-\nu^{\sigma}.
\]
Then we apply the Bernoulli formula
\begin{align*}
\sum_{\mu=1}^{\nu}\mu^{\sigma}&=\frac{\nu^{\sigma+1}}{\sigma+1}+
\frac{B_1}{1}\binom{\sigma}{0}\nu^{\sigma}+
\frac{B_{2}}{2}\binom{\sigma}{1}\nu^{\sigma-1}+\cdots\\
&\quad +\frac{B_{\sigma-1}}{\sigma-1}\binom{\sigma}{\sigma-2}\nu^2+B_{\sigma}\nu,
\end{align*}
where $B_1=1/2$ and $B_{2n+1}=0$ for $n>0$.
\end{proof}
\begin{corollary} \label{cor1}
\[
N_{\sigma,j}(\nu)=\gamma^{\sigma,j}_1\nu^{\sigma+1}+\gamma^{\sigma,j}_{3}\nu^{\sigma-1}+\cdots+
\gamma^{\sigma,j}_{3+2k}\nu^{\sigma-1-2k}+\cdots,\; k=0,1,2,\ldots.
\]
where the last term is proportional to $\nu$ or $\nu^2$ depending on whether $\sigma$ is even or odd.
When $\sigma$ is even, the last term is $2B_{\sigma}\nu$.
\end{corollary}
\begin{proof}
Using (\ref{f45}) and Theorem~\ref{th3} we derive the representation
\[
N_{\sigma,j}(\nu)=\sum_{n=0}^{j}\binom{j}{n}(-1)^n2
\left\{\frac{\nu^{\sigma+1}}{\sigma-n+1}+\frac{1}{2}B_2\binom{\sigma-n}{1}\nu^{\sigma-1}+
\cdots\right\}.
\]
The last term appears in the sum only when $n=0$.
\end{proof}
Since $\sigma$ is even, all monomials in the polynomial $N_{\sigma,j}(\nu)$ have odd degree.
Applying Theorem~\ref{th1} to (\ref{f42}) when $\sigma+1\leq 2s-1$ we obtain
$J_{\sigma,j}=B_{\sigma}$ for $\sigma=2,4,\ldots,2s-2$.

The case $\sigma=2s$ is not covered by Theorem~\ref{th1} because $\sigma+1>2s-1$. It is
treated by formula (\ref{f37}) for $i=j=s$:
\[
\xi_{ij}+\xi_{ji}=2\sum_{\nu=1}^s\alpha_{\nu}\sum_{\mu=0}^{\nu-1}\mu^s(\mu-\nu)^s.
\]

\subsection{Construction of a diagonal positive definite matrix $H$ and skew-symmetric matrix $B_2$}

The diagonal matrix satisfying (\ref{f35}) is positi\-ve definite only if $1\leq s\leq 4$.
The corresponding SBP operator has interior accuracy of order $2s$ and is formally accurate of order $s$
near boundaries. But the overall accuracy is of order $s+1$; see e.g. \cite{26}.

Let $X=Q\begin{bmatrix}\hat{R}\\0\end{bmatrix}$ be a QR factorization of $X$ in (\ref{f29}).
Then the skew-symmetric matrix $B_2$ have the structure $B_2=Q\begin{bmatrix}\Theta_{11}&\Theta_{12}\\\Theta_{21}&\Theta_{22}\end{bmatrix}Q^T$, where
$\begin{bmatrix}\Theta_{11}\\\Theta_{21}\end{bmatrix}=Q^TG\hat{R}^{-1}$, $\Theta_{12}=-\Theta_{21}^T$ and $\Theta_{22}$ is an arbitrary skew-symmetric matrix.
The skew-symmetric matrix $\Theta_{22}$ is of size $(s-1)\times(s-1)$ and therefore depends on $(s-2)(s-1)/2$ free parameters.

The following MATLAB function requires the Symbolic Math Toolbox. The output parameters
\texttt{alpha,h,B} denote the vector $\alpha$ in Theorem~\ref{th1},
the diagonal of $H$ and the matrix $B$, respectively.
The computed matrix $B$ has $\Theta_{22}=0$; see line 14 of the function.
\begin{verbatim}
function [alpha,h,B] = SBP(s)
s = sym(s); r = 2*s;
alpha = -(-1).^(1:s).*cumprod((s:-1:1)./(s+1:2*s))./(1:s);
Cs = toeplitz(alpha(s:-1:1), ...
     [alpha(s) sym(zeros(1,double(s-1)))]);
C = [sym(zeros(double(r-s))); Cs];
X = (r:-1:1)'.^(0:r-1); Y = (0:-1:1-s)'.^(0:s);
F = -C*Y; F(1,:) = F(1,:) + r.^(0:s)/2;
Hankel = X(:,1:s+1)'*F + F'*X(:,1:s+1);
h = Hankel([2:s+1, (s+1)*s+2:(s+1)^2])./(1:r)/X;
G = F - diag(h)*X(:,[1 1:s])*diag(0:s);
[Q,R] = qr(X(:,1:s+1));
T1 = Q'*G/R(1:s+1,:);
T = [T1 [-T1(s+2:r,:)'; sym(zeros(double(r-s-1)))]];
B = Q*T*Q'; B(1,1) = -1/2;
\end{verbatim}
Here is the output for $s=1,2$:
\begin{description}
\item[$s=1$]\quad $\alpha_1=1/2$, $H=\mbox{diag}[1/2,1]$,
$B=\begin{bmatrix}-1/2,&1/2\\-1/2,&0\end{bmatrix}$;
\item[$s=2$]\quad $\alpha_1=2/3$, $\alpha_2=-1/12$, $H=\mbox{diag}[17, 59, 43, 49]/18$,\newline
\[
B=\begin{bmatrix} -1/2,& 59/96, & -1/12, & -1/32\\ -59/96, & 0, & 59/96, &  0\\
  1/12, & -59/96, & 0, & 59/96\\ 1/32, & 0, & -59/96, & 0\end{bmatrix}.
\]
\end{description}

\subsection{Mirror classical SBP operator}

We build an SBP operator that approximates the first derivative on the semi-infinite grid
$x_{\nu}=-\nu h$, $\nu=0,1,\ldots$, in the form of a block-partitioned semi-infinite matrix
\begin{equation} \label{f46}
h^{-1}\begin{bmatrix} \hat{D} & \hat{C} \\ -\hat{H}^{-1}\hat{C}^T & \hat{H}^{-1}\hat{B}\end{bmatrix},
\end{equation}
where
\begin{align*}
\hat{C}&=\begin{bmatrix}\vdots&\vdots\\0&0\\C_s&C_0^T\end{bmatrix}, \quad
C_s=\begin{bmatrix}\alpha_s& 0 &\cdots & \cdots & 0\\
\alpha_{s-1} & \alpha_s & 0 & \cdots & 0\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\vdots& & \ddots& \ddots & 0 \\
\alpha_1 & \cdots & \cdots & \alpha_{s-1} & \alpha_s \end{bmatrix},\\
\hat{D}&=\begin{bmatrix}
\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\vdots\\
0&-\alpha_s & \cdots & -\alpha_1 & 0 & \alpha_1 &\cdots & \alpha_s&0\\
&\ddots&\ddots&-\alpha_2&-\alpha_1 & 0 & \alpha_{1} & \cdots & \alpha_s\\
&&\ddots&\ddots&\vdots&\ddots&\ddots&\ddots&\vdots\\
&&&\ddots&-\alpha_s & \cdots & -\alpha_1 & 0 & \alpha_1\\
&&&&0& -\alpha_s & \cdots & -\alpha_1 & 0
\end{bmatrix}.
\end{align*}
The matrix $C_0$ is the zero $(r-s)\times s$ matrix. The coefficients $\alpha_{\nu}$ in $C_s$ and
in the semi-infinite matrix $\hat{D}$ are defined in Theorem~\ref{th1}.
The semi-infinite block $\hat{D}$ provides $2s$-order accurate approximation of the first derivative
at the inner grid points $x_{\nu}$, $\nu\le -r$.

Let us introduce the $r\times r$ permutation matrix
\[
\Pi=\begin{bmatrix}
&&&1\\
&&\iddots\\
&1\\
1
\end{bmatrix}
\]
which flips a vector $v$ such that $(\Pi v)_i=r+1-i$, $i=1,2,\ldots,r$.
We assume that $r\geq s$ and choose the $r\times r$ matrix $\hat{B}$ from (\ref{f46}) in the form
\begin{equation} \label{f47}
\hat{B}=-\Pi B_1\Pi+\hat{B}_2,\quad B_1=\mbox{diag}(-1/2,0,\ldots,0),\quad \hat{B}_2^T=-\hat{B}_2.
\end{equation}

The algebraic system for approximation of order $\tau$ near the boundary is
\begin{equation} \label{f48}
\hat{H}^{-1}\hat{B}\Pi\mathbf{x}^j-\hat{H}^{-1}\begin{bmatrix}C_s^T\\C_0\end{bmatrix}
\hat{\mathbf{y}}^j=j\Pi\mathbf{x}^{j-1}, \quad j=0,1,\ldots,\tau,
\end{equation}
where $\mathbf{x}=\begin{bmatrix}r&r-1&\ldots&1\end{bmatrix}^T$ and
$\hat{\mathbf{y}}=\begin{bmatrix}1-s&\ldots&-1&0\end{bmatrix}^T$ is the flipped
vector $\mathbf{y}=\begin{bmatrix}0&-1&\ldots&1-s\end{bmatrix}^T$.
Multiplying (\ref{f48}) by $\Pi\hat{H}$ from the left gives the system
\[
\Pi\hat{B}_2\Pi\mathbf{x}^j-B_1\mathbf{x}^j-\Pi
\begin{bmatrix}C_s^T\\C_0\end{bmatrix}\hat{\mathbf{y}}^j=
j\Pi\hat{H}\Pi\mathbf{x}^{j-1}, \quad j=0,1,\ldots,\tau.
\]
As a result, we arrive at the following equations for $j=0,1,\ldots,\tau$
\begin{equation} \label{f49}
\Pi\hat{B}_2\Pi\mathbf{x}^j=-\mathbf{g}_j, \mbox{ where }\;
\mathbf{g}_j=-B_1\mathbf{x}^j-\begin{bmatrix}C_0\\C_s\end{bmatrix}\mathbf{y}^j
-j\Pi\hat{H}\Pi\mathbf{x}^{j-1}.
\end{equation}
Comparing (\ref{f49}) with (\ref{f28}) gives the desired formulas
\begin{equation} \label{f50}
\hat{B}_2=-\Pi B_2\Pi,\quad \hat{H}=\Pi H\Pi.
\end{equation}

\section{Symmetric hyperbolic systems with dissipative boundary conditions}
\label{sec3}

Consider a symmetric hyperbolic system \cite{7,9}
\begin{equation} \label{f51}
A(t,x,y)\frac{\partial}{\partial t}u+B(t,x,y)\frac{\partial}{\partial x}u+
C(t,x,y)\frac{\partial}{\partial y}u+D(t,x,y)u=f(t,x,y)
\end{equation}
for $0\leq t\leq T$, $a\leq x\leq b$, $c\leq y\leq d$. The matrices $A$, $B$ and $C$ are Hermitian that is $A=A^*$, $B=B^*$, $C=C^*$. Moreover, the matrix $A$ is positive definite.
The matrices $A$, $B$, $C$, $D$ and the vector functions $u(t,x,y)$ and $f(t,x,y)$ may have complex entries. The real part of a complex number $z$ is denoted by $\mathsf{Re\,}z$.

\begin{theorem} \label{th4}
Solution of (\ref{f51}) satisfies the matrix identity
\begin{multline} \label{f52}
\frac{\partial}{\partial t}(Au,u)+\frac{\partial}{\partial x}(Bu,u)
 +\frac{\partial}{\partial y}(Cu,u)+2\mathsf{Re}(Du,u)\\=((A_t+B_x+C_y)u,u)
 +2\mathsf{Re}(f,u),
\end{multline}
where $(u,v)=u^*v$ is the Euclidean inner product of vectors $u$ and $v$.
\end{theorem}
The proof of Theorem~\ref{th4} is direct.

Let us introduce the $L^2$ inner product of vector functions $u(t,x,y)$ and $v(t,x,y)$
\[
\langle u,v\rangle=\int_{a}^b\int_c^du^*vdxdy
\]
and the $L^2$ inner products on the boundaries:
\begin{align*}
\langle u,v\rangle|_{x=a}&=\int_c^d(u^*v)|_{x=a}dy,\quad
\langle u,v\rangle|_{x=b}=\int_c^d(u^*v)|_{x=b}dy,\\
\langle u,v\rangle|_{x=a}^{x=b}&=\langle u,v\rangle|_{x=b}-\langle u,v\rangle|_{x=a}.
\end{align*}
\begin{theorem} \label{th5}
Solution of (\ref{f51}) satisfies the energy identity
\begin{multline}\label{f53}
\frac{d}{dt}\langle Au,u\rangle+\langle Bu,u\rangle|_{x=a}^{x=b}
 +\langle Cu,u\rangle|_{y=c}^{y=d}+2\mathsf{Re}\langle Du,u\rangle\\ =\langle(A_t+B_x+C_y)u,u\rangle+2\mathsf{Re}\langle f,u\rangle.
\end{multline}
\end{theorem}
The proof of Theorem~\ref{th5} is direct.

\subsection{Boundary conditions}

Assume that boundary conditions imposed on the solution of the symmetric hyperbolic system (\ref{f51})
are of the form
\begin{equation} \label{f54}
\begin{split}
Z_a(t,y)u&=g_a(t,y)\mbox{ for } x=a,\qquad Z_b(t,y)u=g_b(t,y)\mbox{ for } x=b,\\
Z_c(t,x)u&=g_c(t,x)\mbox{ for } y=c,\qquad Z_d(t,x)u=g_d(t,x)\mbox{ for } y=d,
\end{split}
\end{equation}
where the matrices $Z_a(t,y)$, $Z_b(t,y)$, $Z_c(t,x)$ and $Z_d(t,x)$ have full row rank.
We also assume that $B(t,x,y)$ is nonsingular when $x=a$ or $x=b$, and $C(t,x,y)$ is nonsingular
when $y=c$ or $x=d$. The theory of hyperbolic equations requires the following necessary conditions
on the boundaries:
\begin{description}
\item[$x=a$] the number of rows in $Z_a(t,y)$ is equal to the number of positive eigenvalues
of the matrix $B(t,a,y)$;
\item[$x=b$] the number of rows in $Z_b(t,y)$ is equal to the number of negative eigenvalues
of the matrix $B(t,a,y)$;
\item[$y=c$] the number of rows in $Z_c(t,y)$ is equal to the number of positive eigenvalues
of the matrix $B(t,x,c)$;
\item[$y=d$] the number of rows in $Z_d(t,y)$ is equal to the number of negative eigenvalues
of the matrix $B(t,x,d)$;
\end{description}

The boundary conditions (\ref{f54}) are called \emph{dissipative} if
\begin{equation} \label{f55}
\begin{split}
(B(t,a,y)v,v) & \leq0 \mbox{ whenever } Z_a(t,y)v=0,\\
(B(t,b,y)v,v) & \geq0 \mbox{ whenever } Z_b(t,y)v=0,\\
(C(t,x,c)v,v) & \leq0 \mbox{ whenever } Z_c(t,y)v=0,\\
(C(t,x,d)v,v) & \geq0 \mbox{ whenever } Z_d(t,y)v=0.
\end{split}
\end{equation}
Let the boundary conditions (\ref{f54}) be homogeneous i.e. $g_a\equiv0$, $g_b\equiv0$, $g_c\equiv0$,
$g_d\equiv0$. When the homogeneous boundary conditions are dissipative the following energy inequality
holds:
\begin{equation} \label{f56}
\frac{d}{dt}\langle Au,u\rangle \leq \langle(A_t+B_x+C_y)u,u\rangle-
2\mathsf{Re}\langle Du,u\rangle+2\mathsf{Re}\langle f,u\rangle.
\end{equation}
A symmetric hyperbolic system (\ref{f51}) with dissipative boundary conditions (\ref{f54}), (\ref{f55})
has a unique solution; see e.g. \cite{18,9,2}.

\section{Tensor operations}
\label{sec4}

We are given integer indices $i,j,k$ in the range $1\leq i\leq s_1$, $0\leq j\leq s_2$, $0\leq k\leq s_3$.
The vector $V_i$ is a tensor of order 1, the matrix $M_{ij}$ is a tensor of order 2. We are mainly interested in tensors of order 3 given by the values $U=U_{ijk}$, which are generally complex.
The indices $j$ and $k$ indicate the grid nodes along the $x$ and $y$ axes, respectively.
The index $i$ is usually used for multiplication of $U$ by a matrix coefficient.

\subsection{Vectorization} A tensor $U$ of order 3 can be unfolded into a vector $\vec{U}=\mbox{vec}(U)$ by using the lexicographic ordering
\begin{equation} \label{f57}
\vec{U}_{i+js_1+ks_1(s_2+1)}=U_{ijk}.
\end{equation}
Tensors of order 2 are vectorized similarly, i.e. $\vec{U}_{i+js_1}=U_{ij}$.

\subsection{Inner product}
The inner product of tensors $U$ and $V$ of order 3 is the number
\begin{equation} \label{f58}
(U,V)=(\vec{U},\vec{V})=\sum_{ijk}\bar{U}_{ijk}V_{ijk}
\end{equation}
where the bar over $U_{ijk}$ denotes complex conjugation.

\subsection{Multiplication by a matrix in a certain index mode}
The pro\-duct of the tensor $U_{ijk}$ with a $s_1'\times s_1$ matrix $A$ in mode 1,
that is, with respect to index $i$, is the tensor $V=A\times_1U$ with entries
\begin{equation} \label{f59}
V_{i'jk}=\sum_{i}A_{i'i}U_{ijk}.
\end{equation}
The product of a tensor with a matrix in modes 2 or 3 is defined analogously.

\subsection{Kronecker product of matrices}
The Kronecker product of an $m\times n$ matrix $A$ and a matrix $B$ is the block matrix
\[
A\otimes B=\begin{bmatrix}a_{11}B&a_{12}B&\cdots&a_{1n}B\\a_{21}B&a_{22}B&\cdots&a_{2n}B\\
\vdots&\vdots&\ddots&\vdots\\a_{m1}B&a_{m2}B&\cdots&a_{mn}B\end{bmatrix}.
\]
The Kronecker product satisfies the identity $(A\otimes B)\otimes C=A\otimes(B\otimes C)$.
The Kronecker product of three matrices $A\otimes B\otimes C$ is $(A\otimes B)\otimes C$
by definition. Moreover, $A_1A_2\otimes B_1B_2=(A_1\otimes B_1)(A_2\otimes B_2)$ and
\[
A_1A_2\otimes B_1B_2\otimes C_1C_2=(A_1\otimes B_1\otimes C_1)(A_2\otimes B_2\otimes C_2).
\]
The vectorized equivalent for the operation $V=A\times_1U$ is
\[
\vec{V}=(I\otimes I\otimes A)\vec{U},
\]
where $I$ denotes the identity matrix of suitable size.
Similar expressions hold for multiplication along modes 2 or 3.

\subsection{Multiplication by variable matrix coefficients}
In what follows, we will often operate with sets of matrices $A(t,j,k)$
that depend on a continu\-ous parameter $t$ and integer indices $j,k$.
Then the product $V=AU$ is defined as
\begin{equation} \label{f60}
V_{i'jk}(t)=\sum_i(A(t,j,k))_{i'i}U_{ijk}(t).
\end{equation}
Note that this type of multiplication is always performed in the first dimen\-sion, $i$,
of the tensor $U$.

To rewrite (\ref{f60}) in a vectorized form, we form a block-diagonal matrix
\begin{equation} \label{f61}
\mathcal{A}(t)=\mbox{blockdiag}\{A(t,j,k)\colon 0\leq j\leq s_2,0\leq k\leq s_3\},
\end{equation}
where the diagonal blocks $A(t,j,k)$ are ordered in ascending order with respect to
the linear index $j+k(s_2+1)$. Then the product (\ref{f60}) is given by the vectorized
expression $\vec{V}(t)=\mathcal{A}(t)\vec{U}(t)$.

\subsection{Boundary slices}
Given tensor values $U_{ijk}$ on a rectangular grid $(x_j,y_k)$, the boundary slices of $U$ are
the tensors $U^{[x_0]}$, $U^{[x_{\rm end}]}$, $U^{[y_0]}$, $U^{[y_{\rm end}]}$ of order 2
defined by the entries
\[
U^{[x_0]}_{ik}=U_{i,0,k},\quad U^{[x_{\rm end}]}_{ik}=U_{i,s_2,k},\quad
U^{[y_0]}_{ij}=U_{i,j,0},\quad U^{[y_{\rm end}]}_{ij}=U_{i,j,s_3}.
\]

For a set of matrix coefficients $A(t,j,k)$, we also introduce the boundary subsets
$A^{[x_0]}$, $A^{[x_{\rm end}]}$, $A^{[y_0]}$, $A^{[y_{\rm end}]}$ such that
\begin{align*}
A^{[x_0]}(t,k)&=A(t,0,k),\quad A^{[x_{\rm end}]}(t,k)=A(t,s_2,k),\\
A^{[y_0]}(t,j)&=A(t,j,0),\quad\; A^{[y_{\rm end}]}(t,j)=A(t,j,s_3).
\end{align*}

\section{Energy estimates for semi-discretized symmetric hyperbolic systems}

The main advantage of SBP operators is that they allow one to replicate with little effort the derivation of the energy identities for partial differential equations when performing semi-discrete approximation.

We perform semi-discretization of (\ref{f51}) on the 2-dimensional grid $(x_j,y_k)$ composed from
the equidistant 1-dimensional grids,
\[
a=x_0<x_1<\cdots < x_{s_x}=b\quad \mbox{ and }\quad c=y_0<y_1<\cdots < y_{s_y}=d.
\]
The tensor components $U_{ijk}(t)$ are approximations to the values $u_i(t,x_j,y_k)$.

Recall that the SBP differentiation matrices that approximate the first derivative on a 1-dimensional grid
have the form $\mathsf{D}=P^{-1}Q$, where $P$ is a positive definite diagonal matrix and $Q$ satisfies the SBP property
\[
Q+Q^T=\mathsf{E}=\mbox{diag}[-1,0,0,\ldots,0,1].
\]
We have two SBP operators, $\mathsf{D}_x=P_x^{-1}Q_x$ and $\mathsf{D}_y=P_y^{-1}Q_y$.
For brevity, we will denote
\[
\mathbb{D}_xV=P^{-1}_xQ_x\times_2V\quad \mbox{ and }\quad \mathbb{D}_yV=P^{-1}_yQ_y\times_3V.
\]
The vectorized form of $V(t)=P_x^{-1}Q_x\times_2U(t)$ is
\[
\vec{V}(t)=(I\otimes P_x^{-1}Q_x\otimes I)\vec{U}(t).
\]
Similarly, the vectorized form of $V(t)=P_y^{-1}Q_y\times_2U(t)$ is given by
\[
\vec{V}(t)=(P_y^{-1}Q_y\otimes I\otimes I)\vec{U}(t).
\]

\begin{remark} \label{rem3}
When the matrix coefficient $B$ varies with respect to $x$, the discretization of the term
$B(t,x,y)\frac{\partial}{\partial x}u$ from (\ref{f51}) in the form $B\mathbb{D}_xU$ does not
lead to a discrete analog of the energy identity.
\end{remark}

A working trick is to use the equality $Bu_x=\frac{1}{2}[(Bu)_x+Bu_x-B_xu]$ and discretize
$\frac{1}{2}[(Bu)_x+Bu_x-B_xu]$ instead of $Bu_x$.
Thus, the SBP semi-discretization is applied to the equivalent system
\begin{multline} \label{f62}
A(t,x,y)\frac{\partial}{\partial t}u+\frac{1}{2}\frac{\partial}{\partial x} \left[B(t,x,y)u\right]
+\frac{1}{2}B(t,x,y)\frac{\partial}{\partial x}u-\frac{1}{2}\frac{\partial}{\partial x}B(t,x,y)u\\
\hspace{4.2em}+\frac{1}{2}\frac{\partial}{\partial y} \left[C(t,x,y)u\right]
+\frac{1}{2}C(t,x,y)\frac{\partial}{\partial y}u-\frac{1}{2}\frac{\partial}{\partial y}C(t,x,y)u\\
+D(t,x,y)u=f(t,x,y)
\end{multline}
instead of (\ref{f51}). We do not discretize the derivatives $A_t$, $B_x$ and $C_y$ and assume
for brevity that these matrices are given.

The system (\ref{f62}) is semi-discretized as follows:
\begin{multline} \label{f63}
A(t,j,k)U_t+\frac{1}{2}\mathbb{D}_x\left[B(t,j,k)U\right]
+\frac{1}{2}B(t,j,k)\mathbb{D}_xU-\frac{1}{2}B_x(t,j,k)U\\
+\frac{1}{2}\mathbb{D}_y\left[C(t,j,k)U\right]
+\frac{1}{2}C(t,j,k)\mathbb{D}_yU-\frac{1}{2}C_y(t,j,k)U+D(t,j,k)U\\
=F(t,j,k)+\mathbb{SAT}.
\end{multline}
Dissipative boundary conditions are discretized by the SAT terms in $\mathbb{SAT}$.
The exact expression for $\mathbb{SAT}$ is given later.

Let us introduce the block-diagonal matrices
\begin{equation} \label{f64}
\begin{split}
\mathcal{A} & =\mbox{blockdiag}(A(t,j,k)),\quad \mathcal{A}_t=\mbox{blockdiag}(A_t(t,j,k)),\\
\mathcal{B} & =\mbox{blockdiag}(B(t,j,k)),\quad \mathcal{B}_x=\mbox{blockdiag}(B_x(t,j,k)),\\
\mathcal{C} & =\mbox{blockdiag}(C(t,j,k)),\quad \mathcal{C}_y=\mbox{blockdiag}(C_y(t,j,k)),\\
\mathcal{D} & =\mbox{blockdiag}(D(t,j,k)),\quad \mathcal{P}=P_y\otimes P_x\otimes I.
\end{split}
\end{equation}
The matrices $\mathcal{A}$, $\mathcal{B}$, $\mathcal{C}$, $\mathcal{A}_t$, $\mathcal{B}_x$,
$\mathcal{C}_y$ are hermitian. The matrix $\mathcal{P}$ is diagonal.

\begin{lemma} \label{lem1}
Owing to the special diagonal structure of $\mathcal{P}$ and block structure of matrices in (\ref{f64}),
\begin{align*}
\mathcal{P}\mathcal{A}&=\mathcal{A}\mathcal{P},\quad \mathcal{P}\mathcal{B}=\mathcal{B}\mathcal{P},\quad
\mathcal{P}\mathcal{C}=\mathcal{C}\mathcal{P},\quad \mathcal{P}\mathcal{D}=\mathcal{D}\mathcal{P},\\
\mathcal{P}\mathcal{A}_t&=\mathcal{A}_t\mathcal{P},\quad \mathcal{P}\mathcal{B}_x=\mathcal{B}_x\mathcal{P},\quad
\mathcal{P}\mathcal{C}_y=\mathcal{C}_y\mathcal{P}.
\end{align*}
The commutativity also holds for the diagonal matrices $P_y\otimes \mathsf{E}_x\otimes I$ and
$\mathsf{E}_y\otimes P_x\otimes I$ instead of $\mathcal{P}$.
\end{lemma}
\begin{proof}
Consider, for example, the identity $\mathcal{P}\mathcal{B}=\mathcal{B}\mathcal{P}$. The block with indices $j,k$ in both parts of the identity is the product of $B(t,j,k)$ with the scalar
$P_x(j,j)P_y(k,k)$, where $P_x(j,j)$ and $P_y(k,k)$ are the corresponding diagonal elements in the diagonal matrices $P_x$ and $P_y$. The remaining identities are proved similarly.
\end{proof}

For convenience, we transform (\ref{f63}) equivalently into vectorized form
\begin{multline} \label{f65}
\mathcal{A}\vec{U}_t+\frac{1}{2}(I\otimes \mathsf{D}_x\otimes I)\mathcal{B}\vec{U} +\frac{1}{2}\mathcal{B}(I\otimes \mathsf{D}_x\otimes I)\vec{U}-\frac{1}{2}\mathcal{B}_x\vec{U}\\
+\frac{1}{2}(\mathsf{D}_y\otimes I\otimes I)\mathcal{C}\vec{U} +\frac{1}{2}\mathcal{C}(\mathsf{D}_y\otimes I\otimes I) \vec{U}-\frac{1}{2}\mathcal{C}_y\vec{U}+\mathcal{D}\vec{U} \hspace{3.3em}\\
=\vec{F}+\vec{\mathbb{SAT}}.
\end{multline}

Let us introduce the inner products on the boundaries $x=a$ and $x=b$:
\begin{align}
(\mathcal{B}\vec{U},\vec{U})_{P_y}^{[x_0]}&=\sum_k P_y(k,k)
\sum_{i,i';j=0}B_{ii'}(t,j,k)\bar{U}_{i'jk}U_{ijk} \label{f66}\\
(\mathcal{B}\vec{U},\vec{U})_{P_y}^{[x_{\rm end}]}&=\sum_k P_y(k,k)
\sum_{i,i';j=s_x}B_{ii'}(t,j,k)\bar{U}_{i'jk}U_{ijk} \label{f67}
\end{align}
The inner products on the other boundaries, $(\mathcal{C}\vec{U},\vec{U})_{P_x}^{[y_0]}$ and
$(\mathcal{C}\vec{U},\vec{U})_{P_x}^{[y_{\rm end}]}$, are defined accordingly.
\begin{theorem} \label{th6}
The following energy identity holds for a solution of (\ref{f65}):
\begin{multline} \label{f68}
\frac{d}{dt}(\mathcal{A}\vec{U},\vec{U})_{\mathcal{P}}+
(\mathcal{B}\vec{U},\vec{U})_{P_y}^{[x_{\rm end}]} -(\mathcal{B}\vec{U},\vec{U})_{P_y}^{[x_0]}
+(\mathcal{C}\vec{U},\vec{U})_{P_x}^{[y_{\rm end}]}-(\mathcal{C}\vec{U},\vec{U})_{P_x}^{[y_0]}\\
=([\mathcal{A}_t+\mathcal{B}_x+\mathcal{C}_y]\vec{U},\vec{U})_{\mathcal{P}}-
2\mathsf{Re}(\mathcal{D}\vec{U},\vec{U})_{\mathcal{P}}+
2\mathsf{Re}(\vec{F},\vec{U})_{\mathcal{P}}\\
+(\vec{\mathbb{SAT}},\vec{U})_{\mathcal{P}}+(\vec{U},\vec{\mathbb{SAT}})_{\mathcal{P}}.
\end{multline}
\end{theorem}
\begin{proof}
Introduce the auxiliary matrices
\begin{align*}
\Phi&=\frac{1}{2}(I\otimes \mathsf{D}_x\otimes I)\mathcal{B} +\frac{1}{2}\mathcal{B}(I\otimes \mathsf{D}_x\otimes I) -\frac{1}{2}\mathcal{B}_x\\
\Psi&=\frac{1}{2}(\mathsf{D}_y\otimes I\otimes I)\mathcal{C} +\frac{1}{2}\mathcal{C}(\mathsf{D}_y\otimes I\otimes I) -\frac{1}{2}\mathcal{C}_y
\end{align*}
and rewrite the system in a compact form as
$\mathcal{A}\vec{U}_t+\Phi\vec{U}+\Psi\vec{U}+\mathcal{D}\vec{U}=\vec{F}+\vec{\mathbb{SAT}}$.
Then
\begin{align*}
\frac{d}{dt}(\mathcal{A}\vec{U},\vec{U})_{\mathcal{P}}&=(\mathcal{A}\vec{U}_t,\vec{U})_{\mathcal{P}}+
(\mathcal{A}\vec{U},\vec{U}_t)_{\mathcal{P}}+(\mathcal{A}_t\vec{U},\vec{U})_{\mathcal{P}}\\
&=(-\Phi\vec{U}-\Psi\vec{U}-\mathcal{D}\vec{U}+\vec{F}+\vec{\mathbb{SAT}},\vec{U})_{\mathcal{P}}\\
&\quad +(\vec{U},-\Phi\vec{U}-\Psi\vec{U}-\mathcal{D}\vec{U}+\vec{F}+\vec{\mathbb{SAT}})_{\mathcal{P}}\\
&=-((\mathcal{P}\Phi+\Phi^*\mathcal{P})\vec{U},\vec{U})-((\mathcal{P}\Psi+
\Psi^*\mathcal{P})\vec{U},\vec{U}))\\
&\quad +(\mathcal{A}_t\vec{U},\vec{U})_{\mathcal{P}}-
2\mathsf{Re}(\mathcal{D}\vec{U},\vec{U})_{\mathcal{P}}+
2\mathsf{Re}(\vec{F},\vec{U})_{\mathcal{P}}\\
&\quad +(\vec{\mathbb{SAT}},\vec{U})_{\mathcal{P}}+(\vec{U},\vec{\mathbb{SAT}})_{\mathcal{P}}.
\end{align*}

Using Lemma~\ref{lem1} we derive the equalities
\begin{align*}
\mathcal{P}\Phi+\Phi^T\mathcal{P}&=\frac{1}{2}(P_y\otimes Q_x\otimes I)\mathcal{B}+
\frac{1}{2}\mathcal{B}(P_y\otimes Q_x\otimes I)-\frac{1}{2}\mathcal{P}\mathcal{B}_x\\
&\quad +\frac{1}{2}(P_y\otimes Q_x^T\otimes I)\mathcal{B}+
\frac{1}{2}\mathcal{B}(P_y\otimes Q_x^T\otimes I)-\frac{1}{2}\mathcal{B}_x\mathcal{P}\\
&=\frac{1}{2}(P_y\otimes (Q_x+Q_x^T)\otimes I)\mathcal{B}+
\frac{1}{2}\mathcal{B}(P_y\otimes (Q_x+Q_x^T)\otimes I)-\mathcal{P}\mathcal{B}_x.
\end{align*}
Since $Q_x+Q^T_x=\mathsf{E}_x$ is diagonal, we obtain
\begin{equation} \label{f69}
\mathcal{P}\Phi+\Phi^T\mathcal{P}=
(P_y\otimes \mathsf{E}_x\otimes I)\mathcal{B}-\mathcal{P}\mathcal{B}_x.
\end{equation}
Analogously, $Q_y+Q^T_y=\mathsf{E}_y$ is diagonal and
\begin{equation} \label{f70}
\mathcal{P}\Psi+\Psi^T\mathcal{P}=
(\mathsf{E}_y\otimes P_x\otimes I)\mathcal{C}-\mathcal{P}\mathcal{C}_y.
\end{equation}
As a result,
\begin{multline*}
\frac{d}{dt}(\mathcal{A}\vec{U},\vec{U})_{\mathcal{P}}+
((P_y\otimes \mathsf{E}_x\otimes I)\mathcal{B}\vec{U},\vec{U})+
((\mathsf{E}_y\otimes P_x\otimes I)\mathcal{C}\vec{U},\vec{U})=\\
([\mathcal{A}_t+\mathcal{B}_x+\mathcal{C}_y]\vec{U},\vec{U})_{\mathcal{P}}
-2\mathsf{Re}(\mathcal{D}\vec{U},\vec{U})_{\mathcal{P}}+
2\mathsf{Re}(\vec{F},\vec{U})_{\mathcal{P}}\\
+(\vec{\mathbb{SAT}},\vec{U})_{\mathcal{P}}+(\vec{U},\vec{\mathbb{SAT}})_{\mathcal{P}}.
\end{multline*}

Since $\mathsf{E}_x=\mbox{diag}[-1,0,\ldots,1]$, we obtain the equality
\begin{align*}
(P_y\otimes \mathsf{E}_x\otimes I)\mathcal{B}&=-\mbox{blockdiag}[P_y(k,k)B(j,k)\delta_j^0]\\
&\quad +\mbox{blockdiag}[P_y(k,k)B(j,k)\delta_j^{s_x}],
\end{align*}
where $\delta_j^n=1$ if $j=n$ and $\delta_j^n=0$ if $j\ne n$. Therefore
\[
((P_y\otimes \mathsf{E}_x\otimes I)\mathcal{B}\vec{U},\vec{U})=-(\mathcal{B}\vec{U},\vec{U})_{P_y}^{[x_0]}
+(\mathcal{B}\vec{U},\vec{U})_{P_y}^{[x_{\rm end}]}.
\]
Similarly,
\[
((\mathsf{E}_y\otimes P_x\otimes I)\mathcal{C}\vec{U},\vec{U})=
-(\mathcal{C}\vec{U},\vec{U})_{P_x}^{[y_0]}
+(\mathcal{C}\vec{U},\vec{U})_{P_x}^{[x_{\rm end}]}.
\]
\end{proof}

\subsection{SAT terms on the boundaries}

Let us construct the following penalty terms corresponding to the boundary conditions (\ref{f54}):
\begin{equation} \label{f71}
\begin{split}
\mathsf{B}_a&=-\frac{\sigma_a}{2}Z_a^*\Lambda_a(Z_aU^{[x_0]}-G_a),\\
\mathsf{B}_b&=-\frac{\sigma_b}{2} Z_b^*\Lambda_b(Z_bU^{[x_{\rm end}]}-G_b),\\
\mathsf{B}_c&=-\frac{\sigma_c}{2} Z_c^*\Lambda_c(Z_cU^{[y_0]}-G_c),\\
\mathsf{B}_d&=-\frac{\sigma_d}{2} Z_d^*\Lambda_d(Z_dU^{[y_{\rm end}]}-G_d).
\end{split}
\end{equation}
These functions depend on positive constants $\sigma_a$, $\sigma_b$, $\sigma_c$, $\sigma_d$ and
on positive definite hermitian matrices $\Lambda_a$, $\Lambda_b$, $\Lambda_c$, $\Lambda_d$.
The matrices $\Lambda_a$, $\Lambda_b$, $\Lambda_c$, $\Lambda_d$ are chosen freely based convenience.
In this paper, we choose
\[
\Lambda_a=(Z_aZ_a^*)^{-1},\quad \Lambda_b=(Z_bZ_b^*)^{-1},\quad \Lambda_c=(Z_cZ_c^*)^{-1},\quad
\Lambda_d=(Z_dZ_d^*)^{-1}.
\]

The functions $\mathsf{B}_a$, $\mathsf{B}_b$, $\mathsf{B}_c$, $\mathsf{B}_d$
are extended by zero from the boundaries to the whole domain:
\begin{align*}
&F_a(t,j,k)=\begin{cases}\mathsf{B}_a(t,k),& j=0,\\0,& j>0,\end{cases}
\qquad
F_b(t,j,k)=\begin{cases}\mathsf{B}_b(t,k),& j=s_2,\\0,& j<s_2,\end{cases}\\
&F_a(t,j,k)=\begin{cases}\mathsf{B}_c(t,j),& k=0,\\0,& k>0,\end{cases}
\qquad
F_b(t,j,k)=\begin{cases}\mathsf{B}_d(t,j),& k=s_3,\\0,& j<s_3.\end{cases}
\end{align*}
Such an extension is done solely to comply with formula (\ref{f63}).

The SAT terms take the form
\begin{align}
\mathbb{SAT}_a=P_x^{-1}\times_2 F_a(t,j,k),\quad \mathbb{SAT}_b=P_x^{-1}\times_2F_b(t,j,k), \label{f72}\\
\mathbb{SAT}_c=P_y^{-1}\times_3F_c(t,j,k),\quad \mathbb{SAT}_d=P_y^{-1}\times_3F_d(t,j,k). \label{f73}
\end{align}
Boundary conditions (\ref{f54}) are enforced in (\ref{f63}) by a penalty using the
combi\-nation of SAT terms:
\begin{equation} \label{f74}
\mathbb{SAT}=\mathbb{SAT}_a+\mathbb{SAT}_b+\mathbb{SAT}_c+\mathbb{SAT}_d.
\end{equation}
Since the matrix $P_y$ is diagonal, we use Lemma~\ref{lem1} and derive the identities
\begin{align*}
(\mathbb{SAT}_a,U)_{\mathcal{P}}+(U,\mathbb{SAT}_a)_{\mathcal{P}}=
-\sigma_a(Z_a^*\Lambda_a Z_aU^{[x_0]},U^{[x_0]})_{P_y}\\
+\frac{\sigma_a}{2}(Z_a^*\Lambda_aG_a,U^{[x_0]})_{P_y}+
\frac{\sigma_a}{2}(U^{[x_0]},Z_a^*\Lambda_aG_a)_{P_y}.
\end{align*}
Other SAT terms are treated similarly.

We arrive at the following theorem.
\begin{theorem} \label{th7}
The following energy identity holds for a solution of (\ref{f63}):
\begin{equation} \label{f75}
\begin{split}
\frac{d}{dt}(\mathcal{A}\vec{U},\vec{U})_{\mathcal{P}}&
+(S^{[x_0]}U^{[x_0]},U^{[x_0]})_{P_y}
+(S^{[x_{\rm end}]}U^{[x_{\rm end}]},U^{[x_{\rm end}]})_{P_y}\\
&+(S^{[y_0]}U^{[y_0]},U^{[y_0]})_{P_x}
+(S^{[y_{\rm end}]}U^{[y_{\rm end}]},U^{[y_{\rm end}]})_{P_x}\\
=\;&([\mathcal{A}_t+\mathcal{B}_x+\mathcal{C}_y]\vec{U},\vec{U})_{\mathcal{P}}
-2\mathsf{Re}(\mathcal{D}\vec{U},\vec{U})_{\mathcal{P}}
+2\mbox{\sf Re}(\vec{F},\vec{U})_{\mathcal{P}}\\
&+\sigma_a\mbox{\sf Re}(Z_a^*\Lambda_aG_a,U^{[x_0]})_{P_y}+
\sigma_b\mbox{\sf Re}(Z_b^*\Lambda_bG_b,U^{[x_{\rm end}]})_{P_y}\\
&+\sigma_c\mbox{\sf Re}(Z_c^*\Lambda_cG_c,U^{[y_0]})_{P_x}+
\sigma_d\mbox{\sf Re}(Z_d^*\Lambda_dG_d,U^{[y_{\rm end}}])_{P_x},
\end{split}
\end{equation}
where we use the hermitian matrices
\begin{align}
S^{[x_0]}(t,k)&=-B^{[x_0]}+\sigma_a Z_a^*\Lambda_aZ_a, \;
S^{[x_{\rm end}]}(t,k)=B^{[x_{\rm end}]}+\sigma_b Z_b^*\Lambda_bZ_b, \label{f76}\\
S^{[y_0]}(t,j)&=-B^{[y_0]}+\sigma_c Z_c^*\Lambda_cZ_c, \;\;
S^{[y_{\rm end}]}(t,j)=B^{[y_{\rm end}]}+\sigma_d Z_d^*\Lambda_dZ_d. \label{f77}
\end{align}
\end{theorem}

We will use \emph{strictly dissipative} boundary conditions instead of non-strictly dissipative
conditions as in (\ref{f55}). Strict dissipation assumes the inequalities
\begin{align*}
(-B^{[x_0]}v,v)&>0 \mbox{ for all nonzero vectors }v \mbox{ such that }Z_av=0,\\
(B^{[x_{\rm end}]}v,v)&>0 \mbox{ for all nonzero vectors }v \mbox{ such that }Z_bv=0,\\
(-B^{[y_0]}v,v)&>0 \mbox{ for all nonzero vectors }v \mbox{ such that }Z_cv=0,\\
(B^{[y_{\rm end}]}v,v)&>0 \mbox{ for all nonzero vectors }v \mbox{ such that }Z_dv=0,
\end{align*}

One can choose the parameters $\sigma_a$, $\sigma_b$, $\sigma_c$ and $\sigma_d$
in (\ref{f71}) such that all matrices $S^{[x_0]}$, $S^{[x_{\rm end}]}$, $S^{[y_0]}$, $S^{[y_{\rm end}]}$ in (\ref{f76}) and (\ref{f77}) are non-negative definite.

The SAT method for a weak enforcement of boundary conditions in the method of lines was proposed in \cite{3}.

\subsection{How to choose $\sigma_a$, $\sigma_b$, $\sigma_c$ and $\sigma_d$}

The matrices in (\ref{f76}) and (\ref{f77}) have the general form
\[
S=-\hat{B}+\sigma Z^*\Lambda Z,
\]
where $\hat{B}$ denotes $B^{[x_0]}$, $-B^{[x_{\rm end}]}$, $B^{[y_0]}$, or
$-B^{[y_{\rm end}]}$. Strict dissipation for a boundary condition $ZU=G$
means that $(\hat{B}v,v)<0$ for all nonzero vectors $v$ satisfying
the equation $Zv=0$.

Let $Z^*=Q\hat{R}$ be a QR factorization, where $Q$ is unitary,
$\hat{R}=\begin{bmatrix}R\\0\end{bmatrix}$ and $R$ is nonsingular.
Using the matrix $Q$ we introduce the following block matrices:
\[
Q^*\hat{B}Q=\begin{bmatrix}B_{11}&B_{12}\\B_{21}&B_{22}\end{bmatrix},\quad
Q^*Z^*\Lambda ZQ=\begin{bmatrix}R\Lambda R^*&0\\0&0\end{bmatrix}.
\]
The block $B_{22}$ is negative definite, that is $B_{22}<0$, due to the strict dissipative
boundary condition $ZU=G$. Since $\Lambda=(ZZ^*)^{-1}=(R^*R)^{-1}$,
we have $R\Lambda R^*=I$. As a result,
\begin{equation} \label{f78}
S=-\hat{B}+\sigma Z^*\Lambda Z=
Q\begin{bmatrix}\sigma I-B_{11}&-B_{12}\\-B_{21}&-B_{22}\end{bmatrix}Q^*.
\end{equation}
Since the matrix $-B_{22}$ is positive definite we require non-negativity
of the Schur complement $\sigma I-B_{11}+B_{12}B_{22}^{-1}B_{21}$ to $-B_{22}$.
The Schur complement is non-negative definite if
\begin{equation} \label{f79}
\sigma\ge \lambda_{\max}(B_{11}-B_{12}B_{22}^{-1}B_{21}),
\end{equation}
where $\lambda_{\max}$ denotes the largest eigenvalue of a matrix.

The inequality (\ref{f79}) guarantees that the matrix $S$ is non-negative definite.

\section{SAT terms for the interfaces between spatial domains}

Consider now two adjacent rectangular domains $\Omega=[a,b]\times[c,d]$ and
$\hat{\Omega}=[\hat{a},\hat{b}]\times[c,d]$, where $b=\hat{a}$. The common boundary
consists of the points $(x,y)$ such that $x=b=\hat{a}$, $y\in[c,d]$.
In the continuous case, we set the following interface boundary conditions:
\begin{align*}
u(t,b,y)&=\hat{u}(t,\hat{a},y)\mbox{ in domain } \Omega,\\
\hat{u}(t,\hat{a},y)&=u(t,b,y)\mbox{ in domain } \hat{\Omega}.
\end{align*}

Assume that semi-discretization of a hyperbolic symmetric system with SBP operators
is done separately in $\Omega$ and $\hat{\Omega}$. We also assume dissipative boundary conditions on the external boundaries of $\Omega$ and $\hat{\Omega}$ and discretize these conditions with appropriate SAT terms.
Solution of the semi-discretized problem in $\Omega$ is denoted by $U$ and in $\hat{\Omega}$ by $\hat{U}$.

Following \cite{4}, we construct penalty terms at the interface in the form
\begin{equation} \label{f80}
\mathsf{B}_b(t,k)=\frac{1}{2}B(U^{[x_{\rm end}]}-\hat{U}^{[\hat{x}_0}]),\quad
\hat{\mathsf{B}}_{\hat{a}}(t,k)=-\frac{1}{2}B(\hat{U}^{[\hat{x}_0]}-U^{[x_{\rm end}]}),
\end{equation}
where the matrix $B$ is taken at the corresponding interface nodes.
The penalty functions are extended by zero from the interface boundary to the interior grid points.
The extensions $F_b(t,j,k)$ and $\hat{F}_{\hat{a}}(t,\hat{j},k)$ allow us to
form the corresponding SAT terms
\begin{equation} \label{f81}
\mathbb{SAT}_b=P_x^{-1}\times_2F_b(t,j,k),\qquad
\mathbb{SAT}_{\hat{a}}=P_{\hat{x}}^{-1}\times_2\hat{F}_{\hat{a}}(t,\hat{j},k).
\end{equation}

We add the energy identities provided by Theorem~\ref{th7} separately for each of the two domains:
\begin{equation} \label{f82}
\frac{d}{dt}[(\mathcal{A}\vec{U},\vec{U})_{\mathcal{P}}+
(\hat{\mathcal{A}}\vec{\hat{U}},\vec{\hat{U}})_{\hat{\mathcal{P}}})]+
(BU,U)_{P_y}^{[x_{\rm end}]}
-(B\hat{U},\hat{U})_{P_y}^{[\hat{x}_0]}+\ldots = \Upsilon+\ldots.
\end{equation}
The contribution of $\mathbb{SAT}_b$ and $\mathbb{SAT}_{\hat{a}}$ to the right-hand side of (\ref{f82}) is
\begin{align*}
\Upsilon&=(\mathbb{SAT}_b,U)_{\mathcal{P}}+(\mathbb{SAT}_{\hat{a}},\hat{U})_{\mathcal{\hat{P}}}+
(U,\mathbb{SAT}_b)_{\mathcal{P}}+(\hat{U},\mathbb{SAT}_{\hat{a}})_{\mathcal{\hat{P}}}\\
&=\frac{1}{2}(B(U^{[x_{\rm end}]}-\hat{U}^{[\hat{x}_0]}),U^{[x_{\rm end}]})_{P_y}
+\frac{1}{2}(BU^{[x_{\rm end}]},U^{[x_{\rm end}]}-\hat{U}^{[\hat{x}_0]})_{P_y}\\
&\quad -\frac{1}{2}(B(\hat{U}^{[\hat{x}_0]}-U^{[x_{\rm end}]}),\hat{U}^{[\hat{x}_0]})_{P_y}
-\frac{1}{2}(B\hat{U}^{[\hat{x}_0]},\hat{U}^{[\hat{x}_0]}-U^{[x_{\rm end}]})_{P_y}\\
&=(BU^{[x_{\rm end}]},U^{[x_{\rm end}]})_{P_y}
-(B\hat{U}^{[\hat{x}_0]},\hat{U}^{[\hat{x}_0]})_{P_y}.
\end{align*}
Thus, the boundary terms
$(BU^{[x_{\rm end}]},U^{[x_{\rm end}]})_{P_y}-(B\hat{U}^{[\hat{x}_0]},\hat{U}^{[\hat{x}_0]})_{P_y}$
are comple\-tely compensated by the SAT terms on the interface.

\section{Method of lines with explicit Runge-Kutta methods}

To integrate the semi-discrete system (\ref{f63}) in time, we use explicit Runge-Kutta (RK) methods.
Recall that an $s$-stage explicit Runge-Kutta method is applied to the system of ordinary
differential equations $\frac{dy}{dt}=f(t,y)$, and the integration step from $t_n$ to $t_{n+1}=t_n+\tau$
consists of the following $s$ stages:
\begin{align*}
k_1&=f(t_n,y_n),\\
k_i&=f\left(t_n+c_i\tau,y_n+\tau\sum_{j=1}^{i-1}a_{ij}k_j\right),\quad i=2,\ldots s,\\
y_{n+1}&=y_n+\tau\sum_{i=1}^sb_ik_i.
\end{align*}
The coefficients of the Runge-Kutta methods are conveniently represented by the \emph{Butcher tableau}
\[
\begin{array}{c|c}
c&A\\\hline\\[-2.5ex] & b^T
\end{array}
=
\begin{array}{c|cccc}
c_1&a_{11}&a_{12}&\ldots&a_{1s}\\
c_2&a_{21}&a_{22}&\ldots&a_{2s}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
c_s&a_{s1}&a_{s2}&\ldots&a_{ss}\\
\hline
   &b_{1}&b_{2}&\ldots&b_{s}
\end{array}
\]
Note that the entries $a_{ij}$ for $j\ge i$ in the tableau for explicit Runge-Kutta methods
are zero and are usually not displayed.
A survey of \emph{low-storage} implementations of explicit Runge-Kutta methods is given in \cite{13}.

If an explicit Runge-Kutta method with $s$ stages is applied to a linear system
$\dot{y}=\mathcal{L}y+b$ with constant matrix $\mathcal{L}$ and constant vector $b$,
then it is equivalent to the one-step iteration $y_{n+1}=\mathscr{P}_s(\tau\mathcal{L})y_n+\check{\mathscr{P}}_{s-1}(\tau\mathcal{L})\tau b$,
where $\mathscr{P}_s(z)=\sum_{i=0}^sp_iz^i$ is a polynomial with real coefficients, called the stability function in \cite[Definition 2.1]{11}. For a matrix $A$ and a vector $b$ from the Butcher tableau,
\cite[Proposition 3.1]{11} yields the formula
\[
\mathcal{P}_s(z)=1+z b^T(I-zA)^{-1}\mathbf{1},
\]
where $\mathbf{1}=\begin{bmatrix}1&1&\ldots&1\end{bmatrix}^T$ is a 1-vector.
The polynomial $\check{\mathscr{P}}_{s-1}(z)$ is related with $\mathscr{P}_s(z)$ by the formula $\check{\mathscr{P}}_{s-1}(z)=(\mathscr{P}_s(z)-1)/z$.

An explicit Runge-Kutta method with $s$ stages is accurate of order $r$ if
\begin{equation} \label{f83}
\mathscr{P}_s(z)=\sum_{i=0}^r\frac{z^i}{i!}+\sum_{i=r+1}^sp_iz^i.
\end{equation}
All 3-stage explicit Runge-Kutta methods of order 3 have the polynomials
\[
\mathscr{P}_3(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6},\quad
\check{\mathscr{P}}_2(z)=1+\frac{z}{2}+\frac{z^2}{6}.
\]
All 4-stage explicit Runge-Kutta methods of order 4 have the polynomials
\[
\mathscr{P}_4(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24},\quad
\check{\mathscr{P}}_3(z)=1+\frac{z}{2}+\frac{z^2}{6}+\frac{z^3}{24}.
\]

\subsection{Semi-bounded dynamic systems}

The semi-discretized hyperbolic systems depend on step sizes $h$ in the spatial domain. We therefore pay particular attention to the question of stability, which is uniform even for arbitrarily small $h$.

Following \cite{27,28}, a dynamic system $\dot{u}=\mathcal{L}u$, which depends on parameters $h$, is called \emph{semi-bounded} if there exist a constant $\eta$ independent of $h$ and a positive-definite hermitian matrix $\mathbb{H}$ such that
\[
\mathbb{H}\mathcal{L}+\mathcal{L}^*\mathbb{H}\leq 2\eta\mathbb{H}.
\]
Solution of $\dot{u}=\mathcal{L}u$, subject to arbitrary initial data $u(0)=u_0$, satisfies the estimate
\[
\|u(t)\|_{\mathbb{H}}\leq e^{\eta t}\|u_0\|_{\mathbb{H}},
\]
where $(u,v)_{\mathbb{H}}=u^*\mathbb{H}v$ denotes the inner product with the weight matrix $\mathbb{H}$ and $\|u\|_{\mathbb{H}}=\sqrt{(u,u)_{\mathbb{H}}}$.
Indeed, $\frac{d}{dt}(u,u)_{\mathbb{H}}=(\mathcal{L}u,u)_{\mathbb{H}}+(u,\mathcal{L}u)_{\mathbb{H}}= ([\mathbb{H}\mathcal{L}+\mathcal{L}^*\mathbb{H}]u,u)\leq 2\eta(u,u)_{\mathbb{H}}$ and
$\|u(t)\|_{\mathbb{H}}^2\leq e^{2\eta t}\|u_0\|_{\mathbb{H}}^2$.

When
\begin{equation} \label{f84}
\mathbb{H}\mathcal{L}+\mathcal{L}^*\mathbb{H}\leq 0,
\end{equation}
the operator $\mathcal{L}$ is called \emph{semi-dissipative}; see e.g. \cite{1}.
Solution of $\dot{u}=\mathcal{L}u$ governed by a semi-dissipative operator $\mathcal{L}$
satisfy the monotonic estimate
\begin{equation} \label{f85}
\|u(t)\|_{\mathbb{H}}\leq \|u_0\|_{\mathbb{H}}.
\end{equation}

\subsection{Local stability on the imaginary axis}

\begin{definition} \cite{16} \label{def1}
An explicit Runge-Kutta method is called locally stable on the imaginary axis
if there exists a constant $R_s>0$ such that
\begin{equation} \label{f86}
|\mathscr{P}_s(i\sigma)|\leq1,\quad -R_s\leq\sigma\leq R_s.
\end{equation}
\end{definition}
In other words, the region of absolute stability
$\mathcal{A}_s=\{z\in\mathbb{C}\colon |\mathscr{P}_s|\leq1 \}$
of a locally stable Runge-Kutta method on the imaginary axis must contain a nontrivial interval $[-iR_s,iR_s]$. Condition (\ref{f86}) ensures uniform stability of the RK method applied to the scalar semi-dissipative equation $\dot{u}=i\frac{a}{h}u$, where $i=\sqrt{-1}$, $a$ is a non-zero real number and the parameter $h>0$ can be arbitrarily small.

A precise characterization of local stability on the imaginary axis is given in \cite{16}.
Namely, a general $s$-stage RK method satisfies Definition~\ref{def1} if and only if
\[
\begin{cases}
(-1)^{\frac{r+1}{2}}(p_{r+1}-1)<0,& r \mbox{ is odd},\\
(-1)^{\frac{r+2}{2}}[p_{r+2}-(r+2)p_{r+1}+r+1]<0, & r \mbox{ is even}.
\end{cases}
\]
The condition (\ref{f86}) rules out the forward Euler method as well as
the 2-stage Heun method of order 2 (also known as the modified Euler method).

In the particular cases $s=r=3$ and $s=r=4$ the explicit Runge-Kutta method satisfies Definition~\ref{def1}.

\section{Strong stability of the 3-stage third-order Runge-Kutta methods
for time-invariant semi-dissipative systems}

We begin with time-invariant systems $\frac{d}{dt}u=\mathcal{L}u$ and semi-dissipative matrices
$\mathcal{L}$. To overcome difficulties of extending stability analysis
to time-dependent systems and including lower-order terms, we restrict ourselves to
the so called \emph{strong stability} which requires $u_n$ to decrease monotonically
in an appropriate norm as
\[
\|u_{n+1}\|\leq\|u_n\|,\quad n\in\mathbb{N}.
\]

\begin{remark} \label{rem4}
The semi-discretized symmetric hyperbolic system from (\ref{f63}) and the SAT terms, which are defined
in (\ref{f72}) and (\ref{f73}), yield the system
\begin{multline} \label{f87}
U_t=-\frac{1}{2}A^{-1}\left[\mathbb{D}_xBU+B\mathbb{D}_xU+\mathbb{D}_yCU+C\mathbb{D}_yU\right]\\
+\frac{1}{2}A^{-1}\left[ B_xU+C_yU-DU\right] + A^{-1}F+A^{-1}\mathbb{SAT}.
\end{multline}
The system (\ref{f87}) satisfies the energy identity in terms of the weighted inner product $(\mathcal{A}\vec{U},\vec{V})_{\mathcal{P}}$, as in Theorem~\ref{th7}.
This inner product can be inherited into the method of lines via the following notation:
\begin{equation} \label{f88}
\langle u,v\rangle=(\mathcal{A}u,v)_{\mathcal{P}},\quad \|u\|=\sqrt{\langle u,u\rangle}.
\end{equation}
Under the strict dissipative boundary conditions, the operator $\mathcal{L}$ in (\ref{f87})
is semi-bounded in the inner product (\ref{f88}).
\end{remark}

The stability analysis below is valid for arbitrary inner product $\langle\cdot,\cdot\rangle$ and for arbitrary $\mathcal{L}$ which is semi-dissipative in this inner product, i.e.
\begin{equation} \label{f89}
\langle \mathcal{L}v,v\rangle + \langle v,\mathcal{L}v\rangle \leq 0 \mbox{ for all  vectors }v.
\end{equation}

One step of size $\tau$ with an explicit Runge-Kutta method applied to the system $\frac{d}{dt}u=\mathcal{L}u$ computes a vector $u_+$ such that
\begin{equation} \label{f90}
u_+ = \mathscr{P}(L)u,\quad L = \tau\mathcal{L}.
\end{equation}
The polynomial $\mathscr{P}(z)$ in (\ref{f90}) is the stability function of the RK method.
The square of the norm after one step fulfils
\begin{equation} \label{f91}
\|u_+\|^2-\|u\|^2=\|\mathscr{P}(L)u\|^2-\|u\|^2=
\mathsf{Re}\,\langle(\mathscr{P}(L)-I)u,(\mathscr{P}(L)+I)u\rangle.
\end{equation}

The following procedure is proposed in \cite{20} as a generalization of the proof
given by E.~Tadmor and D.~Levermore in \cite{27}.

\begin{procedure} \label{proc1}
In order to estimate the squared norm (``energy'') after one step (\ref{f90}), proceed as follows:
\begin{enumerate}
\item If there is only one term of lowest order in $L$ on either the left or the right hand side of
the scalar product estimating $\|u_+\|^2-\|u\|^2$, remove this lowest order term using the
semi-boundedness of $L$.
\item If the lowest order terms on the left and right hand side are of equal order, estimate the
remaining terms with respect to $\|L\| = \tau\|\mathcal{L}\|$.
\end{enumerate}
\end{procedure}

\begin{theorem}\cite{27,20} \label{th8}
Consider the stability function of any 3-stage Runge-Kutta method of order 3,
\[
\mathscr{P}_3(z) = 1+z+\frac{1}{2}z^2+\frac{1}{6}z^3.
\]
Suppose that $\mathcal{L}$ is a semi-dissipative matrix in some inner product
$\langle\cdot,\cdot\rangle$, i.e.
\begin{equation} \label{f92}
\mathsf{Re}\,\langle \mathcal{L}v,v\rangle\leq0 \mbox{ for all vectors }v.
\end{equation}
Let the vector norm $\|u\|=\langle u,u\rangle^{1/2}$ be defined by the inner product used
in (\ref{f92}), and $\|\mathcal{L}\|$ be the matrix norm induced by this vector norm.
The Runge-Kutta solution of the time-invariant system $\dot{u}=\mathcal{L}u$
after one step of size $\tau$ is $u_+=\mathscr{P}_3(\tau\mathcal{L})u$.
Then $\|u_+\|\leq\|u\|$ for arbitrary $u$ if $\tau\|\mathcal{L}\|\leq1$.
\end{theorem}
\begin{proof}

\begin{lemma} \label{lem2}
Let $\mathscr{P}(L) = \mathscr{P}_3(L) = I+L+\frac{1}{2}L^2+\frac{1}{6}L^3$. Then for all $u$
\begin{multline*}
\langle(\mathscr{P}(L)-I)u,(\mathscr{P}(L)+I)u\rangle\\=2\langle Lv_1,v_1\rangle +
6\langle v_2,Lv_2\rangle +\langle(-\frac{1}{2}I+\frac{1}{6}L)v_3,\frac{1}{6}(I+L)v_3\rangle,
\end{multline*}
where $v_1=(I+\frac{1}{2}L+\frac{1}{6}L^2)u$, $v_2=\frac{1}{6}(I+L)Lu$, $v_3=L^2u$.
\end{lemma}
\begin{proof}
We apply Procedure~\ref{proc1}.
\begin{align*}
\langle(\mathscr{P}(L)-I)u,&(\mathscr{P}(L)+I)u\rangle\\
&=\langle L(I+\frac{1}{2}L+\frac{1}{6}L^2)u,2(I+\frac{1}{2}L+\frac{1}{4}L^2+\frac{1}{12}L^3)u\rangle\\
&=2\langle L(I+\frac{1}{2}L+\frac{1}{6}L^2)u,(I+\frac{1}{2}L+\frac{1}{6}L^2)u\rangle\\
&\quad +\langle L(I+\frac{1}{2}L+\frac{1}{6}L^2)u,(\frac{1}{6}L^2+\frac{1}{6}L^3)u\rangle
\end{align*}
and
\begin{multline*}
\langle L(I+\frac{1}{2}L+\frac{1}{6}L^2)u,(\frac{1}{6}L^2+\frac{1}{6}L^3)u\rangle
=\langle (L+\frac{1}{2}L^2+\frac{1}{6}L^3)u,L\frac{1}{6}(L+L^2)u\rangle\\
= 6\langle \frac{1}{6}(L+L^2)u,L\frac{1}{6}(L+L^2)u\rangle+
\langle (-\frac{1}{2}L^2+\frac{1}{6}L^3)u,\frac{1}{6}L(L+L^2)u\rangle.
\end{multline*}
Finally,
\begin{align*}
\langle (-\frac{1}{2}L^2+\frac{1}{6}L^3)u,\frac{1}{6}L(L+L^2)u\rangle
=\langle (-\frac{1}{2}I+\frac{1}{6}L)L^2u,\frac{1}{6}(I+L)L^2u\rangle.
\end{align*}
\end{proof}
Lemma 2 and (\ref{f91}) imply the equality
\[
\|u_+\|^2-\|u\|^2=2\mbox{\sf Re}\langle Lv_1,v_1\rangle +6\mbox{\sf Re}\langle v_2,Lv_2\rangle +\mbox{\sf Re}\langle(-\frac{1}{2}I+\frac{1}{6}L)v_3,\frac{1}{6}(I+L)v_3\rangle.
\]
Since $\mbox{\sf Re}\langle Lv,v\rangle\leq0$ for all $v$, we have
$2\mbox{\sf Re}\langle Lv_1,v_1\rangle +6\mbox{\sf Re}\langle v_2,Lv_2\rangle\leq 0$.
It follows that
\begin{align*}
\|u_+\|^2-\|u\|^2&\leq \mbox{\sf Re}\langle(-\frac{1}{2}I+\frac{1}{6}L)L^2u,\frac{1}{6}(I+L)L^2u\rangle\\
&\leq\frac{1}{12}\left[-1+\frac{2}{3}\|L\|+\frac{1}{3}\|L\|^2\right]\|L^2u\|^2.
\end{align*}
The polynomial $\frac{1}{3}x^2+\frac{2}{3}x-1=(x-1)(\frac{1}{3}x+1)$ is not positive for
$0\leq x\leq 1$. Hence $\|u_+\|^2-\|u\|^2\leq 0$ when $\|L\|\leq 1$.
\end{proof}

For completeness, we prove an estimate of the contribution from the vector $b$ in
the system $\frac{d}{dt}u=\mathcal{L}u+b$.
\begin{theorem} \label{th9}
Consider the polynomial
\[
\check{\mathscr{P}}_2(z) = 1+\frac{1}{2}z+\frac{1}{6}z^2.
\]
Under the assumptions of Theorem~\ref{th8}, $\|\check{\mathscr{P}}_2(z)u\|\leq\|u\|$
for arbitrary $u$ when $\tau\|\mathcal{L}\|\leq1$.
\end{theorem}
\begin{proof}

\begin{lemma} \label{lem3}
Let $\mathscr{P}(L) = \check{\mathscr{P}}_2(L) = I+\frac{1}{2}L+\frac{1}{6}L^2$. Then for all vectors $u$
\begin{multline*}
\langle(\mathscr{P}(L)-I)u,(\mathscr{P}(L)+I)u\rangle =\langle Lv_1,v_1\rangle +
\frac{1}{12}\langle(I+\frac{1}{3}L)v_2,(-I+L)v_2\rangle,
\end{multline*}
where $v_1=(I+\frac{1}{3}L)u$, $v_2=Lu$.
\end{lemma}
\begin{proof}
By Procedure~\ref{proc1},
\begin{multline*}
\langle(\mathscr{P}(L)-I)u,(\mathscr{P}(L)+I)u\rangle
=\langle(\frac{1}{2}L+\frac{1}{6}L^2)u,(2I+\frac{1}{2}L+\frac{1}{6}L^2)u\rangle\\
=\langle\frac{1}{2}L(I+\frac{1}{3}L)u,2(I+\frac{1}{4}L+\frac{1}{12}L^2)u\rangle\\
=\langle L(I+\frac{1}{3}L)u,(I+\frac{1}{3}L)u\rangle
 +\langle L(I+\frac{1}{3}L)u,(-\frac{1}{12}L+\frac{1}{12}L^2)u\rangle,
\end{multline*}
\begin{align*}
&\langle L(I+\frac{1}{3}L)u,(-\frac{1}{12}L+\frac{1}{12}L^2)u\rangle
=\frac{1}{12}\langle (I+\frac{1}{3}L)Lu,(-I+L)Lu\rangle.
\end{align*}
\end{proof}
Lemma~\ref{lem3} and (\ref{f91}) give the equality
\[
\|\check{\mathcal{P}}_2(L)u\|^2-\|u\|^2=2\mbox{\sf Re}\langle Lv_1,v_1\rangle +\frac{1}{12}\mbox{\sf Re}\langle(I+\frac{1}{3}L)v_2,(-I+L)v_2\rangle.
\]
Since $\mbox{\sf Re}\langle Lv,v\rangle\leq0$ for all $v$, we have
$2\mbox{\sf Re}\langle Lv_1,v_1\rangle \leq 0$.
It follows that
\begin{align*}
\|\check{\mathscr{P}}_2(L)u\|^2-\|u\|^2&\leq \frac{1}{12}\mbox{\sf Re}\langle(I+\frac{1}{3}L)L^2u,(-I+L)L^2u\rangle\\
&\leq\frac{1}{12}\left[-1+\frac{2}{3}\|L\|+\frac{1}{3}\|L\|^2\right]\|L^2u\|^2.
\end{align*}
The polynomial $\frac{1}{3}x^2+\frac{2}{3}x-1$ is not positive for
$0\leq x\leq 1$. Therefore, $\|\check{\mathscr{P}}_2(L)u\|^2-\|u\|^2\leq 0$ when $\|L\|\leq 1$.
\end{proof}
When applying a 3-stage Runge-Kutta method of order 3 to the system
$\frac{d}{dt}y=\mathcal{L}y+b$ with semi-dissipative $\mathcal{L}$, we have the strong
stability estimate
\begin{equation} \label{f93}
\|u_{n+1}\|\leq\|u_n\|+\tau\|b\|,\quad n\in\mathbb{N}.
\end{equation}

\section{Stability of SSPRK3 for time-dependent systems}

When the matrix $\mathcal{L}$ depends on time $t$, the Runge-Kutta scheme can no longer be written as the iteration $y_{n+1}=\mathscr{P}_s(\tau\mathcal{L})y_n$, where $\mathscr{P}_s(z)$ is a truncated series for the exponential. Moreover, for example, the various 3-stage third-order Runge-Kutta schemes are no longer equivalent.

We demonstrate how to investigate stability of time-dependent systems using the SSPRK3 Runge-Kutta method \cite{22} as an example. Other RK methods can be treated similarly; see, e.g., \cite{24}.
The SSPRK3 method is defined by the stages
\begin{equation} \label{f94}
\begin{split}
k_1&=f(t_n,y_n),\quad k_2=f(t_n+\tau,y_n+\tau k_1),\\
k_3&=f(t_n+\frac{\tau}{2},y_n+\frac{\tau}{4}(k_1+k_2)),\\
y_{n+1}&=y_n+\frac{\tau}{6}\left[k_1+k_2+4k_3\right].
\end{split}
\end{equation}
One step of SSPRK3 reduces to the iteration $y_{n+1}=\mathcal{R}(t_n)y_n$ with
\begin{equation} \label{f95}
\begin{split}
\mathcal{R}(t_n)= I&+\frac{\tau}{6}\left[\mathcal{L}(t_n)+\mathcal{L}(t_{n+1})+4\mathcal{L}(t_{n+1/2})\right]\\
&+ \frac{\tau^2}{6}\left[\mathcal{L}(t_{n+1})\mathcal{L}(t_{n})+
   \mathcal{L}(t_{n+1/2})\mathcal{L}(t_{n})+\mathcal{L}(t_{n+1/2})\mathcal{L}(t_{n+1})\right]\\
&+ \frac{\tau^3}{6}\mathcal{L}(t_{n+1/2})\mathcal{L}(t_{n+1})\mathcal{L}(t_{n}).
\end{split}
\end{equation}

\begin{proposition} \label{prop1}
Suppose that $\mathcal{L}$ satisfies the Lipschitz condition, that is,
there exists positive $\zeta$ such that for all $t'$, $t''$
\begin{equation}\label{f96}
\|\mathcal{L}(t')-\mathcal{L}(t''\|\leq\zeta|t'-t''|.
\end{equation}
Then
\begin{align*}
\|\mathcal{R}(t_n)-\mathscr{P}_3(\tau\mathcal{L}(t_n))\|&\leq
\tau\frac{\tau\zeta}{4}\left[1+(1+\tau\|\mathcal{L}(t_n)\|)^2\right]\\
&\quad +\tau^2\frac{(\tau\zeta)^2}{12}(1+\tau\|\mathcal{L}(t_n)\|).
\end{align*}
\end{proposition}
\begin{proof}
If we denote $\delta_1=\mathcal{L}(t_{n+1})-\mathcal{L}(t_n)$ and
$\delta_{1/2}=\mathcal{L}(t_{n+1/2})-\mathcal{L}(t_n)$, then (\ref{f95})
allows us to derive that
\begin{align*}
\mathcal{R}(t_n)&=\mathscr{P}_3(\tau\mathcal{L}(t_n))+
\frac{\tau}{6}\left(\delta_1+4\delta_{1/2}\right)\\
&\quad +\frac{\tau^2}{6}\left[\delta_1\mathcal{L}(t_n)+\delta_{1/2}\mathcal{L}(t_n)+
\delta_{1/2}\mathcal{L}(t_n)+\mathcal{L}(t_n)\delta_1+\delta_{1/2}\delta_1 \right]\\
&\quad +\frac{\tau^3}{6}[\delta_{1/2}\mathcal{L}(t_n)+\mathcal{L}(t_n)\delta_1+
\delta_{1/2}\delta_1]\mathcal{L}(t_n).
\end{align*}
Due to (\ref{f96}), $\|\delta_1\|\leq \tau\zeta$ and $\|\delta_{1/2}\|\leq \frac{\tau}{2}\zeta$.
Therefore,
\begin{align*}
\|\mathcal{R}(t_n)&-\mathscr{P}_3(\tau\mathcal{L}(t_n))\|\\
&\leq\frac{\tau}{6}3\tau\zeta +\frac{\tau^2}{6}3\tau\zeta\|\mathcal{L}(t_n)\|+
\frac{\tau^2}{6}\frac{1}{2}\tau^2\|\mathcal{L}(t_n)\|^2\\
&\quad +\frac{\tau^3}{6}\frac{3}{2}\tau\zeta\eta^2+
\frac{\tau^3}{6}\frac{1}{2}\tau^2\zeta^2\|\mathcal{L}(t_n)\|\\
&=\tau\frac{\tau\zeta}{4}\left[1+(1+\tau\|\mathcal{L}(t_n)\|)^2\right]+
\tau^2\frac{(\tau\zeta)^2}{12}(1+\tau\|\mathcal{L}(t_n)\|).
\end{align*}
\end{proof}
Suppose that $\tau\|\mathcal{L}(t_n)\|\leq1$ and the matrix $\mathcal{L}(t_n)$ is
semi-dissipative. Then $\|\mathscr{P}_3(\tau\mathcal{L}(t_n))\|\leq 1$ by Theorem~\ref{th8}.
By Proposition~\ref{prop1}
\[
\|\mathcal{R}(t_n)\|\leq 1+\tau\frac{5\tau\zeta}{4}+\tau^2\frac{(\tau\zeta)^2}{6}.
\]
The structure of (\ref{f87}) implies that
\[
\zeta=\frac{\zeta_1}{h_x}+\frac{\zeta_2}{h_y}+\zeta_3
\]
with positive constants $\zeta_1$, $\zeta_2$ and $\zeta_3$.
The step-sizes $h_x$ and $h_y$ can be arbitrarily small. Hence
$\tau\zeta=\frac{\tau}{h_x}\zeta_1+\frac{\tau}{h_y}\zeta_2+\tau\zeta_3$. Since
$\tau\|\mathcal{L}(t_n)\|\leq1$, there is a constant $C>0$ such that following inequality holds:
\[
\tau\leq C\min(h_x,h_y).
\]
It follows that $\tau\zeta=C(\zeta_1+\zeta_2)+\tau\zeta_3$ and
\begin{align}
\|\mathcal{R}(t_n)\|&\leq 1+\tau\left\{\frac{5}{4}[C(\zeta_1+\zeta_2)+\tau\zeta_3]
+\frac{1}{6}\tau[C(\zeta_1+\zeta_2)+\tau\zeta_3]^2\right\} \label{f97} \\
&=1+\tau C(\zeta_1+\zeta_2)+O(\tau^2)=1+O(\tau). \label{f98}
\end{align}
The bound (\ref{f98}) ensures Lax-Richtmyer stability \cite{19}.

\subsection{Lower order terms}

The above derivation of bound (\ref{f98}) does not take into account the lower order terms
\[
\mathcal{T}(t)=\frac{1}{2}A^{-1}\left[ B_xU+C_yU-DU\right]
\]
in (\ref{f87}). They are included by adding $\mathcal{T}(t)$ to all matrices $\mathcal{L}(t)$ in (\ref{f95}):
\begin{multline} \label{f99}
\tilde{\mathcal{R}}(t_n)= I+\frac{\tau}{6}\left[\mathcal{L}(t_n)+\mathcal{T}(t_n)+
\mathcal{L}(t_{n+1})+\mathcal{T}(t_{n+1})+\ldots\right]\\
 +\frac{\tau^2}{6}\left[(\mathcal{L}(t_{n+1})+\mathcal{T}(t_{n+1}))
   (\mathcal{L}(t_{n})+\mathcal{T}(t_{n}))+\ldots\right]\\
 +\frac{\tau^3}{6}(\mathcal{L}(t_{n+1/2})+\mathcal{T}(t_{n+1/2}))
(\mathcal{L}(t_{n+1})+\mathcal{T}(t_{n+1}))
   (\mathcal{L}(t_{n})+\mathcal{T}(t_{n})).
\end{multline}
For convenience, we denote
$\theta=\max\{\|\mathcal{L}(t)\|\colon t\in\{t_n,t_{n+1/2},t_{n+1}\}\}$.
Let us suppose that
\[
\tau\theta\leq 1 \mbox{ and } \|\mathcal{T}(t)\|\leq\eta \mbox{ for }
t\in\{t_n,t_{n+1/2},t_{n+1}\}.
\]
Then
\begin{align}
\|\tilde{\mathcal{R}}(t_n)-\mathcal{R}(t_n)\|&\leq\frac{\tau}{6}3\eta+
\frac{\tau^2}{6}3[2\theta\eta+\eta^2]+\frac{\tau^3}{6}3[3\theta^2\eta+3\theta\eta^2+\eta^3] \notag\\
&\leq\frac{\tau}{2}\left[\eta+2\eta+\tau\eta^2+3\eta+3\tau\eta^2+\tau^2\eta^3\right] \notag\\
&= \frac{\tau}{2}\left[6\eta+4\tau\eta^2+\tau^2\eta^3\right]=3\tau\eta+O(\tau^2). \label{f100}
\end{align}
Combining (\ref{f100}) with (\ref{f98}) gives the stability bound
$\|\tilde{\mathcal{R}}(t_n)\|\leq 1+O(\tau)$.

\section{On the stability of the classical Runge-Kutta method}

The classical fourth-order Runge-Kutta method \cite{17} has 4 stages:
\begin{equation} \label{f101}
\begin{split}
k_1&=f(t_n,y_n), \hspace{5.9em} k_2=f(t_n+\frac{\tau}{2},y_n+\frac{\tau}{2}k_1),\\
k_3&=f(t_n+\frac{\tau}{2},y_n+\frac{\tau}{2}k_2),\quad
k_4=f(t_n+\tau,y_n+\tau k_3),\\
y_{n+1}&=y_n+\frac{\tau}{6}\left[k_1+2k_2+2k_3+k_4\right].
\end{split}
\end{equation}

The method has the stability function $\mathscr{P}_4(z)=1+z+\frac{z^2}{2!}+\frac{z^3}{3!}+\frac{z^4}{4!}$.
In \cite{24} it is shown that $\|\mathscr{P}_4(\tau\mathcal{L})\|>1$, in the 2-norm, for the semi-dissipative matrix
\[
\mathcal{L}=-\begin{bmatrix}1 & 2 & 2\\ 0 & 1 & 2\\0& 0& 1 \end{bmatrix}
\]
and for arbitrarily small $\tau>0$.

However, $\|\mathscr{P}_4^2(\tau\mathcal{L})\|\leq1$ for arbitrary semi-dissipative $\mathcal{L}$ and sufficiently small $\tau>0$. This fact was first proved in \cite{24}. Below we give an alternative proof which is based on Procedure~\ref{proc1}. Note that the proof of Lemma~\ref{lem4} is relatively simple when using a computer program for symbolic computation.

\begin{lemma} \label{lem4}
Let $\mathscr{P}(L) = I+L+\frac{1}{2}L^2+\frac{1}{6}L^3+\frac{1}{24}L^4$. Then for all vectors $u$
\begin{align*}
\langle(\mathscr{P}^2(L)-I)u,(\mathscr{P}^2(L)+I)u\rangle&
=\langle Lv_1,v_1\rangle +3\langle v_2,Lv_2\rangle+8\langle Lv_3,v_3\rangle \\
 &\quad +\langle(p_5(L)L^3u,q_5(L)L^3u\rangle,
\end{align*}
where $\langle\cdot,\cdot\rangle$ is any inner product and
\begin{align*}
v_1&=(2I+2L+\frac{4}{3}L^2+\frac{2}{3}L^3+\frac{1}{4}L^4+\frac{5}{72}L^5+
\frac{1}{72}L^6+\frac{1}{576}L^7)u,\\
v_2&=(\frac{2}{3}I+\frac{2}{3}L+\frac{5}{12}L^2+\frac{13}{72}L^3+\frac{1}{18}L^4+\frac{7}{576}L^5+
\frac{1}{576}L^6)Lu,\\ v_3&=(\frac{1}{12}I+\frac{1}{8}L+\frac{1}{12}L^2+\frac{19}{576}L^3+\frac{5}{576}L^4+
\frac{1}{576}L^5)L^2u,\\
p_5(z)&=\frac{1}{12}+\frac{1}{8}z+\frac{1}{12}z^2+\frac{19}{576}z^3+\frac{5}{576}z^4+
\frac{1}{576}z^5,\\
q_5(z)&=-\frac{1}{3}-\frac{1}{4}z-\frac{1}{12}z^2-\frac{1}{72}z^3-\frac{1}{576}z^4+
\frac{1}{576}z^5.
\end{align*}
\end{lemma}
Since the matrix $L$ is semi-dissipative, we obtain the inequality
\[
\mathsf{Re}\langle Lv_1,v_1\rangle +3\langle v_2,Lv_2\rangle+8\langle Lv_3,v_3\rangle\leq0.
\]
Then we have $\langle(p_5(L)v,q_5(L)v\rangle\leq r(\|L\|)\|v\|^2$, where the polynomial $r(z)$
is obtained from the product $p_5(z)q_5(z)$ by replacing its coefficients with their absolute values,
but excluding the coefficient $-1/36$ for zero degree:
\begin{align*}
r(z)&=-\frac{1}{36}+\frac{z}{16}+\frac{19z^2}{288}+\frac{25z^3}{576}+\frac{23z^4}{1152}
+\frac{31z^5}{4608}\\
&\quad +\frac{z^6}{648}+\frac{59z^7}{331776}+\frac{z^8}{55296}+\frac{z^9}{82944}
+\frac{z^{10}}{331776}.
\end{align*}
Since $r(x)<0$ for $x\in[0,0.3147]$, we arrive at the desired stability bound
\begin{equation} \label{f102}
\|\mathscr{P}_4^2(\tau\mathcal{L})\|<1 \mbox{ when } \|\tau\mathcal{L}\|\leq 0.3147.
\end{equation}


\bigskip

\begin{thebibliography}{1}

\bibitem{1}
F. Achleitner, A. Arnold, A. J\"{u}ngel,
\href{https://doi.org/10.1007/978-3-031-65195-3_1}
{\it Necessary and sufficient conditions for strong stability of explicit Runge-Kutta methods},
In Carlen, E., Gon\c{c}alves, P., Soares, A.J. (eds.)
From Particle Systems to Partial Differential Equations. PSPDE 2022,
Springer Proc. in Math. \& Statistics, vol. \textbf{465}, Springer, Cham, (2024), 1--21.

\bibitem{2}
S. Benzoni-Gavage, D. Serre,
\href{https://doi.org/10.1093/acprof:oso/9780199211234.001.0001}
{\it Multidimensional Hyperbolic Partial Differential Equa\-tions.
 First-order Systems and Applications}, Oxford Univ. Press, NY, 2007.

\bibitem{3}
M.H. Carpenter, D. Gottlieb, S. Abarbanel,
\href{https://doi.org/10.1006/jcph.1994.1057}
{\it Time-stable boundary conditions for finite-difference schemes solving
hyperbolic systems: methodology and application to high-order compact schemes},
J. Comput. Phys., \textbf{111} (1994), 220--236.

\bibitem{4}
M.H. Carpenter, J. Nordstr\"{o}m, D. Gottlieb,
\href{https://doi.org/10.1006/jcph.1998.6114}
{\it A stable and conservative interface treatment of arbitrary spatial accuracy},
J. Comput. Phys., \textbf{148} (1999), 341--365.

\bibitem{5}
D.C. Del Rey Fernandez, P.D. Boom, D.W. Zingg,
\href{https://doi.org/10.1016/j.jcp.2014.01.038}
{\it A generalized framework for nodal first derivative summation-by-parts operators},
J. Comput. Phys., \textbf{266} (2014), 214--239.

\bibitem{6}
D.C. Del Rey Fernandez, J.E. Hicken, D.W. Zingg,
\href{https://doi.org/10.1016/j.compfluid.2014.02.016}
{\it Review of summation-by-parts operators with simultaneous approximation terms
for the numerical solution of partial differential equations},
Computers \& Fluids, \textbf{95} (2014), 171--196.

\bibitem{7}
K.O. Friedrichs,
\href{https://doi.org/10.1002/cpa.3160070206}
{\it Symmetric hyperbolic linear differential equations},
Comm. Pure Appl. Math., \textbf{7} (1954), 345--392.

\bibitem{8}
J. Glaubitz, J. Nordstr\"{o}m, P. \"{O}ffner,
\href{https://doi.org/10.1016/j.jcp.2024.112889}
{\it Summation-by-parts operators for general function spaces},
SIAM J. Numer. Anal., \textbf{61}:2 (2023), 733--754.

\bibitem{9}
S.K. Godunov,
{\it Equations of Mathematical Physics}, 2nd ed., Nauka, Moscow, 1979.

\bibitem{10}
B. Gustafsson, H.-O. Kreiss, J. Oliger,
\href{https://doi.org/10.1002/9781118548448}
{\it Time Dependent Problems and Difference Methods}, John Wiley \& Sons, Inc., 1995.

\bibitem{11}
E. Hairer, G. Wanner,
\href{https://doi.org/10.1007/978-3-642-05221-7}
{\it Solving Ordinary Differential Equations II. Stiff and Differen\-tial-Algebraic Problems},
2nd ed., Springer, Berlin, 1996.

\bibitem{12}
R.A. Horn, C.R. Johnson,
\href{https://doi.org/10.1017/CBO9780511810817}
{\it Matrix Analysis}, 2nd ed., Cambridge University Press, NY, 2013.

\bibitem{13}
C.A. Kennedy, M.H. Carpenter, R.M. Lewis,
\href{https://doi.org/10.1016/S0168-9274(99)00141-5}
{\it Low-storage, explicit Runge-Kutta schemes for the compressible Navier-Stokes equations},
Appl. Numer. Math., \textbf{35} (2000), 177--219.

\bibitem{14}
H.-O. Kreiss, G. Scherer,
\href{https://doi.org/10.1016/B978-0-12-208350-1.50012-1}
{\it Finite element and finite difference methods for hyperbolic partial differential equations},
in Mathematical Aspects of Finite Elements in Partial Differential Equations, Academic Press, NY
(1974), 195--212.

\bibitem{15}
H.-O. Kreiss, G. Scherer,
\href{https://user.it.uu.se/~kot25782/fulltext/Scherer77.pdf}
{\it On the existence of energy estimates for difference approximations for hyperbolic systems},
manuscript, 1977.

\bibitem{16}
H.-O. Kreiss, G. Scherer,
\href{https://doi.org/10.1137/0729041}
{\it Method of lines for hyperbolic differential equations},
SIAM J. Numer. Anal., \textbf{29}:3 (1992), 640--646.

\bibitem{17}
M. Kutta,
\href{https://archive.org/details/zeitschriftfrma12runggoog/page/434/mode/2up?view=theater}
{\it Beitrag zur n\"{a}herungsweisen Integration totaler Differentialgleichungen},
Zeitschrift f\"{u}r Mathematik und Physik, \textbf{46} (1901), 435--453.

\bibitem{18}
P.D. Lax, R.S. Phillips,
\href{https://doi.org/10.1002/cpa.3160130307}
{\it Local boundary conditions for dissipative symmetric linear differential operators},
Comm. Pure Appl. Math., \textbf{13} (1960), 427--455.

\bibitem{19}
P.D. Lax, R.D. Richtmyer,
\href{https://doi.org/10.1002/cpa.3160090206}
{\it Survey of the stability of linear finite difference equations},
Comm. Pure Appl. Math., \textbf{9}:2 (1956), 267--293.

\bibitem{20}
H. Ranocha, P. \"{O}ffner,
\href{https://doi.org/10.1007/s10915-017-0595-4}
{\it $L_2$ stability of explicit Runge-Kutta schemes},
J. Sci. Comput., \textbf{75} (2018), 1040--1056.

\bibitem{21}
G. Scherer,
\href{https://user.it.uu.se/~kot25782/fulltext/Scherer77.pdf}
{\it On Energy Estimates for Difference Approximations to Hyperbolic Partial Differential Equations},
PhD Thesis, Uppsala Univ., Uppsala, Sweden, 1977.

\bibitem{22}
Chi-Wang Shu, S. Osher,
\href{https://doi.org/10.1016/0021-9991(88)90177-5}
{\it Efficient implementation of essentially non-oscillatory shock-capturing schemes},
J. Comput. Phys., \textbf{77}, (1988), 439--471.

\bibitem{23}
B. Strand,
\href{https://doi.org/10.1006/jcph.1994.1005}
{\it Summation by parts for finite difference approximations for $d/dx$},
J. Comput. Phys., \textbf{110}:1 (1994), 47--67.

\bibitem{24}
Z. Sun, Chi-Wang~Shu,
\href{https://dx.doi.org/10.4310/AMSA.2017.v2.n2.a3}
{\it Stability of the fourth order Runge-Kutta method for time-dependent partial differential equations},
Annals Math. Sci. Appl., \textbf{2}:2 (2017), 255--284.

\bibitem{25}
M. Sv\"{a}rd, J. Nordstr\"{o}m,
\href{https://doi.org/10.1016/j.jcp.2014.02.031}
{\it Review of summation-by-parts schemes for initial-boundary value problems},
J. Comput. Phys., \textbf{268} (2014), 17--38.

\bibitem{26}
M. Sv\"{a}rd, J. Nordstr\"{o}m,
\href{https://doi.org/10.1016/j.jcp.2019.07.018}
{\it On the convergence rates of energy-stable finite-difference schemes},
J. Comput. Physics, \textbf{397}, (2019), 108819.

\bibitem{27}
E. Tadmor,
\href{https://www.math.umd.edu/~tadmor/pub/linear-stability/Tadmor.SIAM-02.pdf}
{\it From semi-discrete to fully discrete: stability of Runge-Kutta schemes by the energy method. II}
In ``Collected Lectures on the Preservation of Stability under Discretization'', Lecture
Notes from Colorado State Univ. Conf., Fort Collins, CO, 2001 (D. Estep and S. Tavener, eds.)
Proc. Appl. Math. \textbf{109}, SIAM (2002), 25--49.

\bibitem{28}
E. Tadmor,
\href{https://doi.org/10.48550/arXiv.2312.15546}
{\it Runge-Kutta methods are stable},
arXiv:2312.15546v1, 2023.

\end{thebibliography}

\end{document} 