# Stability of finite difference schemes for the hyperbolic initial boundary value problem by winding number computations

Benjamin Boutin, Pierre Le Barbenchon, Nicolas Seguin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Benjamin Boutin, Pierre Le Barbenchon, Nicolas Seguin. Stability of finite difference schemes for the hyperbolic initial boundary value problem by winding number computations. Annales de la Faculté des Sciences de Toulouse. Mathématiques., inPress, 10.48550/arXiv.2302.01410. hal-03970775

HAL Id: hal-03970775

https://hal.science/hal-03970775

Submitted on 2 Feb 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Stability of finite difference schemes for the hyperbolic initial boundary value problem by winding number computations

Benjamin Boutin, Pierre Le Barbenchon, Nicolas Seguin

## ► To cite this version:

Benjamin Boutin, Pierre Le Barbenchon, Nicolas Seguin. Stability of finite difference schemes for the hyperbolic initial boundary value problem by winding number computations. Annales de la Faculté des Sciences de Toulouse. Mathématiques., 2025, 10.5802/afst.1809. hal-03970775v2

HAL Id: hal-03970775

https://hal.science/hal-03970775v2

Submitted on 7 Oct 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# STABILITY OF FINITE DIFFERENCE SCHEMES FOR THE HYPERBOLIC INITIAL BOUNDARY VALUE PROBLEM BY WINDING NUMBER COMPUTATIONS

BENJAMIN BOUTIN, PIERRE LE BARBENCHON, AND NICOLAS SEGUIN

**ABSTRACT.** In this paper, we present a numerical strategy to check the strong stability (or GKS-stability) of one-step explicit finite difference schemes for the one-dimensional advection equation with an inflow boundary condition. The strong stability is studied using the Kreiss-Lopatinskii theory. We introduce a new tool, the intrinsic Kreiss-Lopatinskii determinant, which possesses the same regularity as the vector bundle of discrete stable solutions. By applying standard results of complex analysis to this determinant, we are able to relate the strong stability of numerical schemes to the computation of a winding number, which is robust and cheap. The study is illustrated with the O3 scheme and the fifth-order Lax-Wendroff (LW5) scheme together with a reconstruction procedure at the boundary.

**AMS classification:** 65M12, 65M06, 30E10

**Keywords:** IBVP, Kreiss-Lopatinskii determinant, GKS-stability, finite difference methods, winding number

## 1. INTRODUCTION

1.1. Motivations and assumptions. The purpose of this work is to establish an efficient numerical strategy to determine whether a given finite difference method on the half line is stable or not. We work on an approximation of the rightgoing linear transport equation set on the positive real axis:

$$ \left\{ \begin{array}{ll} \partial_t u + a \partial_x u = 0, & t \ge 0, x \ge 0, \\ u(t, 0) = g(t), & t \ge 0, \\ u(0, x) = f(x), & x \ge 0, \end{array} \right. \qquad (1) $$

where $u(t,x) \in \mathbb{R}$ is the unknown, $f$ an initial datum at time $t = 0$, $g$ is a prescribed physical boundary datum at the point $x = 0$ which corresponds to the inflow boundary because the velocity $a$ is assumed to be positive $a > 0$.

At the discrete level, we consider explicit one-step finite difference methods of the form

$$ U_j^{n+1} = \sum_{k=-r}^{p} a_k U_{j+k}^{n}, \qquad (2) $$

with integers $r, p \ge 1$ and $a_p, a_{-r}$ non zero. The case where $p = 0$ or $r = 0$ will be discussed in Section 1.2. Here, the unknown of the scheme $U_j^n$ is expected to approximate the quantity $u(n\Delta t, j\Delta x)$. The time step $\Delta t > 0$ and the space step $\Delta x > 0$ are usually chosen with respect to some CFL condition $\lambda = a\Delta t/\Delta x \le \lambda_{CFL}$ discussed later on.

Throughout this paper we denote $\mathbb{S} = \{z \in \mathbb{C}, |z| = 1\}$ the unit circle, $\mathbb{D} = \{z \in \mathbb{C}, |z| < 1\}$ the open unit disk, $\mathcal{U} = \{z \in \mathbb{C}, |z| > 1\}$ the associated exterior domain and $\bar{\mathcal{U}} = \{z \in \mathbb{C}, |z| \ge 1\}$ its closure. For $n < m$, the notation $\llbracket n : m \rrbracket$ is for the set $\{k \in \mathbb{N}, n \le k \le m\}$.

<footnote>Date: February 2, 2023.</footnote>

<footnote>This work has been partially supported by ANR project NABUCO, ANR-17-CE40-0025 and by Centre Henri Lebesgue, program ANR-11-LABX-0020-0.</footnote>
<figure><img src="image_2.png" /></figure>

Annales de la Faculté
des Sciences de Toulouse
MATHÉMATIQUES

BENJAMIN BOUTIN, PIERRE LE BARBENCHON AND NICOLAS SEGUIN

*Stability of finite difference schemes for the hyperbolic initial boundary value problem by winding number computations*

Tome XXXIV, n° 1 (2025), p. 193-224.

https://doi.org/10.5802/afst.1809

© les auteurs, 2025.

Les articles des *Annales de la Faculté des Sciences de Toulouse* sont mis à disposition sous la license Creative Commons Attribution (CC-BY) 4.0 http://creativecommons.org/licenses/by/4.0/

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Publication membre du centre
Mersenne pour l'édition scientifique ouverte
http://www.centre-mersenne.org/
e-ISSN : 2258-7519
As a central idea in numerical analysis, the Lax equivalence theorem [25] asserts that a consistent scheme is convergent if and only if it is stable. Therefore, all along the paper, only consistent numerical schemes are considered and the discussion concentrates only on their stability issues. The Cauchy-stability for the space-periodic problem is handled with the Fourier symbolic analysis, the so-called Von-Neumann stability analysis (see [8] and [9]) and makes use of the symbol $\gamma$. The symbol associated with the scheme (2) is defined, for $\xi \in \mathbb{R}$, by

$$ \gamma(\xi) = \sum_{k=-r}^{p} a_k e^{ik\xi}. \qquad (3) $$

**Assumption (H1).** The scheme (2) is Cauchy-stable, meaning that the symbol $\gamma$ satisfies $|\gamma(\xi)| \le 1$ for all $\xi \in \mathbb{R}$.

When dealing with discrete schemes set over the full line $j \in \mathbb{Z}$, the algebraic characterization of the Cauchy-stability follows classically from the Fourier analysis but in the scalar case, it reduces to a geometric property concerning the *symbol curve* $\Gamma$ which is a closed complex parametrized curve defined by

$$ \Gamma = \{\theta \in [0, 2\pi] \mapsto \gamma(\theta)\}. $$

This curve enables a geometric interpretation of the Cauchy-stability assumption (H1) which can be seen as the inclusion $\Gamma \subset \bar{\mathbb{D}}$. The stability condition (H1) can be easily illustrated graphically in the complex plane. In some sense, our goal is to extend this kind of graphical study when including the numerical boundary conditions.

For solving the Initial Boundary Value Problem (IBVP) (1) with the discrete scheme (2), $r$ additional ghost points are needed to take into account the left boundary condition and to fully define the discrete approximation. We assume that the values at these ghost points are obtained from a linear combination of the first values of the solution close to the boundary and at the same time step. More clearly, the considered numerical schemes reads

$$ U_j^{n+1} = \sum_{k=-r}^{p} a_k U_{k+j}^n, \quad j \in \mathbb{N}, n \in \mathbb{N}, \qquad (4) $$

$$ U_j^n = \sum_{k=0}^{m-1} b_{j,k} U_k^n + g_j^n, \quad j \in [-r: -1], n \in \mathbb{N}, \qquad (5) $$

$$ U_j^0 = f_j, \quad j \in \mathbb{N}, \qquad (6) $$

where the integer $m$ satisfies $p+r \le m$, $(f_j)_j$ are approximations of the initial condition $f$ and $(g_j^n)_{n,j}$ are numerical data related to the boundary datum $g$ and possibly its derivatives (see for instance the example in Section 4.3). The assumption $p+r \le m$ is not restrictive since some of the coefficients $b_{j,k}$ are possibly zero.

In order to define the stability on $\ell^2(\mathbb{N})$ and for the sake of convenience in the Kreiss-Lopatinskii determinant formulation (see Definition 10), the explicit use of the $r$ ghost points $U_j^n$, for $j \in [-r: -1]$, can be avoided by substituting the $r$ boundary condition (5) into the recurrence formula (4) for $j \in [0: r-1]$. After straightforward calculations, the boundary part reads also under the form

$$ \mathcal{U}_r^{n+1} = \mathcal{B}\mathcal{U}_m^n + \mathcal{G}^n \qquad (7) $$

where we denote

$$ \mathcal{U}_r^{n+1} = \begin{pmatrix} U_0^{n+1} \\ \vdots \\ U_{r-1}^{n+1} \end{pmatrix}, \quad \mathcal{U}_m^n = \begin{pmatrix} U_0^n \\ \vdots \\ U_{m-1}^n \end{pmatrix}, \quad \mathcal{G}^n = \begin{pmatrix} a_{-r} & \cdots & a_{-1} \\ & \ddots & \\ 0 & & a_{-r} \end{pmatrix} \begin{pmatrix} g_{-r}^n \\ \vdots \\ g_{-1}^n \end{pmatrix} \in \mathcal{M}_r(\mathbb{C}). $$
# Stability of finite difference schemes for the hyperbolic initial boundary value problem by winding number computations (*)

BENJAMIN BOUTIN (1),
PIERRE LE BARBENCHON (2) AND NICOLAS SEGUIN (3)

**ABSTRACT.** — In this paper, we present a numerical strategy to check the strong stability (or GKS-stability) of one-step explicit finite difference schemes for the one-dimensional advection equation with an inflow boundary condition. The strong stability is studied using the Kreiss–Lopatinskii theory. We introduce a new tool, the intrinsic Kreiss–Lopatinskii determinant, which possesses the same regularity as the vector bundle of discrete stable solutions. By applying standard results of complex analysis to this determinant, we are able to relate the strong stability of numerical schemes to the computation of a winding number, which is robust and cheap. The study is illustrated with the O3 scheme and the fifth-order Lax–Wendroff (LW5) scheme together with a reconstruction procedure at the boundary.

**RÉSUMÉ.** — Dans cet article, nous présentons une stratégie numérique permettant de vérifier la stabilité forte (ou stabilité GKS) des schémas de différences finies explicites à un pas pour l'équation d'advection unidimensionnelle avec une condition de bord entrante. La stabilité forte est étudiée à l'aide de la théorie de Kreiss–Lopatinskii. Nous introduisons un nouvel outil, le déterminant intrinsèque de Kreiss–Lopatinskii, qui possède la même régularité que le fibré vectoriel des solutions stables. En lui appliquant les résultats usuels de l'analyse complexe, nous sommes en mesure de ramener l'étude de la stabilité forte au calcul d'indice d'un lacet, procédure fiable et peu coûteuse. L'étude est illustrée avec les schémas O3 et de Lax–Wendroff d'ordre cinq (LW5) avec une condition de bord obtenue par reconstruction.

**Keywords:** IBVP, Kreiss–Lopatinskii determinant, GKS-stability, finite difference methods, winding number.

2020 Mathematics Subject Classification: 65M12, 65M06, 30E10.

This work has been partially supported by ANR project NABUCO, ANR-17-CE40-0025 and by Centre Henri Lebesgue, program ANR-11-LABX-0020-0.

Article proposé par Sylvie Benzoni.

<footnote>(*) Reçu le 10 mai 2023, accepté le 29 novembre 2023.</footnote>

<footnote>(1) Univ Rennes, CNRS, IRMAR - UMR 6625, 35000 Rennes, France — benjamin.boutin@univ-rennes.fr</footnote>

<footnote>(2) Univ Rennes, CNRS, IRMAR - UMR 6625, 35000 Rennes, France — pierre.le-barbenchon@ens-rennes.fr</footnote>

<footnote>(3) Antenne Inria de l'Université de Montpellier, IMAG, CNRS, université de Montpellier, France — nicolas.seguin@inria.fr</footnote>
Here, the matrix $\mathcal{B} \in \mathcal{M}_{r,m}(\mathbb{C})$ encodes the boundary treatment in another way. It corresponds to the boundary part of the quasi-Toeplitz matrix form of the scheme used by Beam and Warming [2]. In the detail, the explicit relationship between $\mathcal{B}$ and $\mathcal{B}$ is as follows:

$$ \mathcal{B} = \begin{pmatrix} a_{-r} & \cdots & a_{-1} \\ \vdots & \ddots & \vdots \\ 0 & \cdots & a_{-r} \end{pmatrix} B + \begin{pmatrix} a_0 & \cdots & a_p & 0 & \cdots & \cdots & 0 \\ \vdots & \ddots & \vdots & \ddots & \ddots & \ddots & \vdots \\ a_{-r+1} & \cdots & a_0 & \cdots & a_p & 0 & \cdots & 0 \end{pmatrix} \in \mathcal{M}_{r,m}(\mathbb{C}) \quad (8) $$

with the notation

$$ B = \begin{pmatrix} b_{-r,0} & \cdots & \cdots & b_{-r,m-1} \\ \vdots & & & \vdots \\ b_{-1,0} & \cdots & \cdots & b_{-1,m-1} \end{pmatrix} \in \mathcal{M}_{r,m}(\mathbb{C}). $$

For example, for the very naive scheme $U_j^{n+1} = \frac{U_{j-1}^n + U_{j+1}^n}{2}$ and the boundary condition $U_{-1}^n = \frac{U_0^n + U_1^n}{2}$, we obtain $B = (\frac{1}{2} \quad \frac{1}{2})$ and $\mathcal{B} = (\frac{1}{4} \quad \frac{3}{4})$.

This class of boundary conditions, (5) or (7), encompasses the Dirichlet and Neumann extrapola- tion procedures, for example, refer to the work of Goldberg [16]. This class also takes into account the more general simplified inverse Lax-Wendroff procedure analyzed by Vilar and Shu [33] in the framework of central compact schemes, and Li, Shu and Zhang for the advection equation [27] and for diffusion equations [28]. We will focus on the so-called reconstruction technique for the bound- ary condition, which enables to deal with a boundary which is not superposed with a grid point (presented by Dakin, Després and Jaouen [10] and also in Section 4.3) in our numerical examples. Other treatments at the boundary exist, as for example absorbing boundary conditions [12] and [11], or transparent boundary conditions [1] and [7], however, in general, they do not enter the present framework.

For finite difference schemes applied to discrete IBVP's, the stability study is a principal issue and is the subject of different approaches. For example, Beam and Warming [2] study the spectral properties of the Toeplitz or quasi-Toeplitz representation of the scheme. In the same spirit, the computation of the spectral radius of the truncated (i.e. finite dimensional) quasi-Toeplitz matrix may provide significant information for the power boundedness of the method. This is the method used for example by Dakin, Després and Jaouen [10]. This strategy is sometimes called eigenvalue spectrum visualization method, especially by Li, Shu and Zhang [26, 27, 28]. In Section 4, we will compare this latter approach with our own strategy presented hereafter for the O3 scheme in the case of reconstruction boundary conditions. Our strategy is based on the so-called GKS-stability theory introduced by Gustafsson, Kreiss and Sundström [19] which handles the discrete IBVP (4)-(5)-(6) with a zero initial data. The reader can refer to the work by Wu [34] and Coulombel [6] for more recent developments on semigroup estimates in order to deduce the stability of the discrete IBVP (4)-(5)-(6) with non zero initial data from the GKS-stability. The notion of GKS-stability (or also called strong stability) for the boundary problem makes use of the following discrete norms:

$$ \|U_j\|_{\Delta t}^2 = \sum_{n=0}^{+\infty} \Delta t |U_j^n|^2 \quad \text{and} \quad \|U\|_{\Delta x, \Delta t}^2 = \sum_{n=0}^{+\infty} \sum_{j=0}^{+\infty} \Delta t \Delta x |U_j^n|^2. $$

The so-called strong stability, or GKS-stability, is defined by:

**Definition 1** (Strong stability). The scheme (4)-(5)-(6) is strongly stable if, for $(f_j) = 0$, there exist $C > 0$ and $\alpha_0$, such that for all $\alpha > \alpha_0$, for all boundary data $(g_j^n)$, for all $\Delta x > 0$, for all $n \in \mathbb{N}$, the solution satisfies

$$ \sum_{j=-r}^{-1} \|e^{-\alpha n \Delta t} U_j\|_{\Delta t}^2 + \left( \frac{\alpha - \alpha_0}{\alpha \Delta t + 1} \right) \|e^{-\alpha n \Delta t} U\|_{\Delta x, \Delta t}^2 \le C \sum_{j=-r}^{-1} \|e^{-\alpha n \Delta t} g_j\|_{\Delta t}^2. \quad (9) $$
# 1. Introduction

## 1.1. Motivations and assumptions

The purpose of this work is to establish an efficient numerical strategy to determine whether a given finite difference method on the half line is stable or not. We work on an approximation of the rightgoing linear transport equation set on the positive real axis:

$$ egin{cases} \partial_t u + a \partial_x u = 0, & t \ge 0, x \ge 0, \\ u(t, 0) = g(t), & t \ge 0, \\ u(0, x) = f(x), & x \ge 0. \end{cases} \tag{1.1} $$

Here $u(t, x) \in \mathbb{R}$ is the unknown, $f$ an initial datum at time $t=0$, and $g$ is a physical boundary datum prescribed at the point $x=0$ which corresponds to the inflow boundary because the velocity $a$ is assumed to be positive $a > 0$.

At the discrete level, we consider explicit one-step finite difference methods of the form

$$ U_j^{n+1} = \sum_{k=-r}^{p} a_k U_{j+k}^n, \tag{1.2} $$

with integers $r, p \ge 1$ and $a_p, a_{-r}$ non zero. The case with $p=0$ or $r=0$ is recalled in Section 1.2. Here, the unknown of the scheme $U_j^n$ is expected to approximate the quantity $u(n\Delta t, j\Delta x)$. The time step $\Delta t > 0$ and the space step $\Delta x > 0$ are usually chosen with respect to some CFL condition $\lambda = a\Delta t/\Delta x \le \lambda_{CFL}$ discussed later on.

**Notation 1.1.** — Throughout this paper we denote $\mathbb{S} = \{z \in \mathbb{C}, |z| = 1\}$ the unit circle, $\mathbb{D} = \{z \in \mathbb{C}, |z| < 1\}$ the open unit disk, $\mathcal{U} = \{z \in \mathbb{C}, |z| > 1\}$ the associated exterior domain and $\bar{\mathcal{U}} = \{z \in \mathbb{C}, |z| \ge 1\}$ its closure. For $n < m$, the notation $\llbracket n : m \rrbracket$ is for the set $\{k \in \mathbb{N}, n \le k \le m\}$.

As a central idea in numerical analysis, the Lax equivalence theorem [30] asserts that a consistent scheme is convergent if and only if it is stable. All of the numerical schemes we consider in this paper are consistent and the discussion only focuses on their stability issues. In addition, we don't address any convergence or accuracy studies. Such results arise easily from usual techniques, once stability is ensured. The Cauchy-stability for the space-periodic problem is handled with the Fourier symbolic analysis, the so-called Von Neumann stability analysis (see [9] and [10]) and makes use of the symbol $\gamma$. The symbol associated with the scheme (1.2) is defined, for $\xi \in \mathbb{R}$, by

$$ \gamma(\xi) = \sum_{k=-r}^{p} a_k e^{ik\xi}. \tag{1.3} $$
We warn the reader that $\|e^{-\alpha n \Delta t} U_j\|_{\Delta t}^2$ is here an abuse of notation to describe $\sum_{n=0}^{+\infty} \Delta t e^{-2\alpha n \Delta t} |U_j^n|^2$ and similarly for $\|e^{-\alpha n \Delta t} U\|_{\Delta x, \Delta t}^2$. This stability definition admits a similar but continuous form for the solutions to continuous hyperbolic PDE's [3]. Namely, it provides some a priori estimates that are useful for a general analysis of such problems.

The following Kreiss theorem [23] expresses a necessary and sufficient condition for the strong stability. We provide hereafter a condensed formulation of this theorem, obtained from [19, Thm 5.1] combined with [18, Lem 13.1.4] or with [17, Def 2.23].

**Theorem 2 (Kreiss).** The following statements are equivalent:

(i) The scheme (4)-(5)-(6) is strongly stable in the sense of Definition 1.

(ii) The Uniform Kreiss-Lopatinskii Condition is satisfied.

The Uniform Kreiss-Lopatinskii Condition corresponds to the absence of zeros for the so-called Kreiss-Lopatinskii determinant $\Delta_{KL}$ that we present here by this informal definition:

$$ \Delta_{KL}(z) = \det(\mathfrak{B}e_1(z), \dots, \mathfrak{B}e_r(z)) \qquad (10) $$

where $(e_1(z), \dots, e_r(z))$ is an explicit basis of the linear space of the $\ell^2(\mathbb{N})$-stable solutions of the $\mathcal{Z}$-transform of the interior equation (4) and $\mathfrak{B}$ is an encoding of the $\mathcal{Z}$-transform of the boundary equation (5). For a proper definition of this determinant, the reader can look at Definition 10 or the book by Gustafsson, Kreiss and Oliger [18]. Before going on, let us provide some comments to a particular case we already studied.

1.2. **The case of totally upwind schemes and summary of [4].** The present article is a non trivial extension of our previous work [4] that deals with the restricted case of totally upwind schemes. Totally upwind schemes are schemes of the form (2) with $p = 0$ if $a > 0$ or $r = 0$ if $a < 0$. Without loss of generality, we restrict here the discussion to the case $p = 0$ since flipping the indices may turn a case to the other. In this section, we summarize the result of [4] and introduce the novelty of the present work. The first step of the analysis conducted in [4] is based on the introduction of the intrinsic Kreiss-Lopatinskii determinant:

$$ \Delta(z) = \frac{\det(\mathfrak{B}e_1(z), \dots, \mathfrak{B}e_r(z))}{\det(e_1(z), \dots, e_r(z))} \qquad (11) $$

using the same informal notation as in (10). Under appropriate assumptions, an explicit formula for the intrinsic Kreiss-Lopatinskii determinant is obtained:

$$ \forall |z| \ge 1, \quad \Delta(z) = (-1)^{r(m-r)} \det C(z) \left( \frac{a-r}{a_0-z} \right)^{m-r} \qquad (12) $$

where $\det C(z)$ is an computable polynomial in $z$ depending only on the coefficients $(a_j)_{j=-r}^0$ and on $\mathfrak{B}$. Thanks to this result, we prove that $\Delta$ is holomorphic on $\bar{U}$. Note that this property may be wrong as long as the standard Kreiss-Lopatinskii determinant is concerned.

Applying the residue theorem to $\Delta$, we develop a numerical strategy to count the number of zeros of the Kreiss-Lopatinskii determinant in $\mathcal{U}$. By Theorem 2 (Kreiss), we conclude that if $\mathrm{Ind}_{\Delta(\mathfrak{S})}(0) < r$ then the scheme is not stable where $\mathrm{Ind}_{\Delta(\mathfrak{S})}(0)$ is the notation for the winding number of 0 with respect to the Kreiss-Lopatinskii curve $\Delta(\mathfrak{S})$. This result allows us to establish an efficient and practical method (see Method 19 of [4] or Method 15 of the present paper) to study the stability of a scheme with boundary. It provides sharp results for the solution $(U_j^n)_{n \in \mathbb{N}} \in \ell^2(\mathbb{N})$ to the problem set on the half line $\mathbb{N}$. In particular and contrary to numerical investigations of stability which are based on the computation of the spectral radius, no arbitrary truncation of (quasi-) Toeplitz matrices is needed. In return, a problem set on a bounded space domain needs, for a whole convergence study, superposition techniques for truncated data, as used in [7] and [5]. This feature restricts mainly the study to explicit scheme.
ASSUMPTION 1.2. — The scheme (1.2) is Cauchy-stable, meaning that the symbol $\gamma$ satisfies $|\gamma(\xi)| \le 1$ for all $\xi \in \mathbb{R}$.

When dealing with discrete schemes set over the full line $j \in \mathbb{Z}$, the algebraic characterization of the Cauchy-stability traditionally follows from the Fourier analysis. In the scalar case, it reduces to a geometric property concerning the symbol curve $\Gamma$ which is a closed complex parametrized curve defined by

$$ \Gamma = \{\theta \in [0, 2\pi] \mapsto \gamma(\theta)\}. $$

This curve enables a geometric interpretation of the Cauchy-stability assumption 1.2 reformulated equivalently using the inclusion $\Gamma \subset \bar{\mathbb{D}}$. The stability condition 1.2 can easily be illustrated in the complex plane. In some sense, our goal is to extend this kind of graphical study when including the numerical boundary conditions.

For solving the Initial Boundary Value Problem (IBVP) (1.1) with the discrete scheme (1.2), $r$ additional ghost points are needed to take into account the left boundary condition and to fully define the discrete approximation. We assume that the values at these ghost points are obtained from a linear combination of the first values of the solution close to the boundary and at the same time step. More clearly, the considered numerical schemes read

$$ U_j^{n+1} = \sum_{k=-r}^{p} a_k U_{k+j}^n, \quad j \in \mathbb{N}, n \in \mathbb{N}, \qquad (1.4) $$

$$ U_j^n = \sum_{k=0}^{m-1} b_{j,k} U_k^n + g_j^n, \quad j \in [-r: -1], n \in \mathbb{N}, \qquad (1.5) $$

$$ U_j^0 = f_j, \quad j \in \mathbb{N}, \qquad (1.6) $$

where the integer $m$ satisfies $p + r \le m$, $(f_j)_j$ approximates the initial condition $f$, and $(g_j^n)_{n,j}$ are numerical data related to the boundary datum $g$ (possibly involving also its derivatives, see for instance the example in Section 4.3). The assumption $p + r \le m$ is not restrictive since some of the coefficients $b_{j,k}$ are possibly zero.

In order to define the stability on $\ell^2(\mathbb{N})$ and for the sake of convenience in the Kreiss–Lopatinskii determinant formulation (see Definition 2.9), the explicit use of the $r$ ghost points $U_j^n$, for $j \in [-r: -1]$, can be avoided by substituting the $r$ boundary condition (1.5) into the recurrence formula (1.4) for $j \in [0: r-1]$. After straightforward calculations, the boundary part reads also under the form

$$ \mathcal{U}_r^{n+1} = \mathcal{B}\mathcal{U}_m^n + \mathcal{G}^n \qquad (1.7) $$
In the present article, we extend to $p \ge 1$ the connection between the winding number of 0 and the stability of the scheme. Indeed, even if we do not have an explicit formula as (12), the holomorphic property of the intrinsic Kreiss-Lopatinskii determinant is sufficient to have the same efficient method to study the stability. This approach, using the winding number, is robust since instead of finding zeros of an algebraic curve, it only requires the computation of a winding number, which is an integer, to count the number of zeros. One can mention the work of Thuné [32] who develops a numerical method to check the GKS-stability. He looks for the precise location of the zeros of the Kreiss-Lopatinskii determinant approximating the roots of some parameterized characteristic polynomial equations which is significantly different with our work.

### 1.3. Outline of the paper.
After constructing the intrinsic Kreiss-Lopatinskii determinant in Section 2 and Section 3, we see that, in such a general case, the lack of an explicit formula for $\Delta(z)$ does not preclude holomorphic properties (see Theorem 12). From there, we obtain the following stability criterion: if $\mathrm{Ind}_{\Delta(S)}(0) < r$ then the scheme is not stable (see Corollary 14). We prove these results by the use of Hermite interpolation and residue theorem. To compute the Kreiss-Lopatinskii determinant numerically, in Section 4, we use an easy-to-use formulation of it which is, in some sense, close to the explicit formulation (12). Moreover, Section 4 gathers the numerical procedure to draw the Kreiss-Lopatinskii curve, several examples, and numerical experiments for illustrating the efficiency of the proposed strategy.

## 2. KREISS-LOPATINSKII DETERMINANTS

In this section, we introduce the Kreiss-Lopatinskii determinant, a usual tool to check the Uniform Kreiss-Lopatinskii Condition. Then we define the intrinsic Kreiss-Lopatinskii determinant, namely a reshaping of the previous one, which is more convenient in practice and has better properties than the classical Kreiss-Lopatinskii determinant: holomorphicity, continuity, independence on the basis...

### 2.1. Stable subspace $\mathcal{E}^s(z)$ and matrix representation.
First, we study the solutions to the interior equation:

$$ U_j^{n+1} = \sum_{k=-r}^{p} a_k U_{k+j}^n, \quad j \in \mathbb{N}, n \in \mathbb{N}. \qquad (13) $$

To study this equation, the $\mathcal{Z}$-transform (see [15, Lesson 40]) is applied. This transformation is defined for $(x_n)_{n \in \mathbb{N}} \in \ell^2(\mathbb{N})$ such that $x_0 = 0$ and $z \in \mathcal{U}$ by $\tilde{x}(z) = \sum_{n \ge 0} z^{-n} x_n$. The previous equation then reads

$$ z \tilde{U}_j(z) = \sum_{k=-r}^{p} a_k \tilde{U}_{j+k}(z), \quad j \in \mathbb{N}, z \in \mathcal{U}. \qquad (14) $$

To solve the linear recurrence equation (14), let us introduce the following characteristic equation where $z$ plays the role of a parameter and $\kappa$ is the indeterminate:

$$ z\kappa^r = \sum_{k=-r}^{p} a_k \kappa^{r+k}. \qquad (15) $$

This equation is nothing but the discrete dispersion relation of the finite difference scheme (13), with frequency parameter $\kappa$ in space and $z$ in time. It is formally obtained by looking for solutions to the interior equation (13) having the form $U_j^n = z^n \kappa^j$.

In the spirit of a classic result by Hersh [22], the following lemma provides a property of separation for the roots with respect to the unit circle.

**Lemma 3 (Hersh).** Assume (H1). For $z$ in the unbounded connected component of $\mathbb{C} \setminus \Gamma$, (1) there is no root of the characteristic equation (15) on $\mathbb{S}$,
where we denote

$$ 
\begin{gathered}
\mathfrak{U}_r^{n+1} = \begin{pmatrix} U_0^{n+1} \\ \vdots \\ U_{r-1}^{n+1} \end{pmatrix}, \quad \mathfrak{U}_m^n = \begin{pmatrix} U_0^n \\ \vdots \\ U_{m-1}^n \end{pmatrix}, \\
\mathcal{G}^n = \begin{pmatrix} a_{-r} & \cdots & a_{-1} \\ \vdots & \ddots & \vdots \\ 0 & a_{-r} \end{pmatrix} \begin{pmatrix} g_{-r}^n \\ \vdots \\ g_{-1}^n \end{pmatrix} \in \mathcal{M}_r(\mathbb{C}).
\end{gathered}
$$

Here, the matrix $\mathcal{B} \in \mathcal{M}_{r,m}(\mathbb{C})$ encodes the boundary treatment in another way. It corresponds to the boundary part of the quasi-Toeplitz matrix form of the scheme used by Beam and Warming [2]. In the detail, the explicit relationship between $\mathcal{B} \in \mathcal{M}_{r,m}(\mathbb{C})$ and $\mathcal{B} \in \mathcal{M}_{r,m}(\mathbb{C})$ is as follows:

$$ \mathcal{B} = \begin{pmatrix} a_{-r} & \cdots & a_{-1} \\ \vdots & \ddots & \vdots \\ 0 & a_{-r} \end{pmatrix} B + \begin{pmatrix} a_0 & \cdots & a_p & 0 & \cdots & \cdots & 0 \\ \vdots & \ddots & \vdots & \vdots & \ddots & \cdots & \vdots \\ a_{-r+1} & \cdots & a_0 & \cdots & a_p & 0 & \cdots & 0 \end{pmatrix} \quad (1.8) $$

with the notation

$$ B = \begin{pmatrix} b_{-r,0} & \cdots & \cdots & b_{-r,m-1} \\ \vdots & & & \vdots \\ b_{-1,0} & \cdots & \cdots & b_{-1,m-1} \end{pmatrix}. $$

The relation (1.8) is invertible since the coefficient $a_{-r}$ is supposed to be non zero. For example, for the very naive scheme $U_j^{n+1} = \frac{U_{j-1}^n + U_{j+1}^n}{2}$ and the boundary condition $U_{-1}^n = \frac{U_0^n + U_1^n}{2}$, we obtain $B = (\frac{1}{2} \quad \frac{1}{2})$ and $\mathcal{B} = (\frac{1}{4} \quad \frac{3}{4})$.

This class of boundary conditions, (1.5) or (1.7), encompasses the Dirichlet and Neumann extrapolation procedures, see, for example, the work of Goldberg [18]. This class also takes into account the more general simplified inverse Lax-Wendroff procedure analyzed by Vilar and Shu [40] in the framework of central compact schemes, and Li, Shu and Zhang for the advection equation [34] and for diffusion equations [35]. We will focus on the so-called reconstruction technique for the boundary condition, which enables to deal with a boundary which is not superposed with a grid point (presented by Dakin, Després and Jaouen [11] and also in Section 4.3) in our numerical examples. Other treatments at the boundary exist, as for example absorbing boundary conditions [13] and [14], or transparent boundary conditions [1] and [8], however, in general, they do not enter the present framework.

For finite difference schemes applied to discrete IBVP's, the stability study is a principal issue and is the subject of different approaches. For example, Beam and Warming [2] study the spectral properties of the Toeplitz
(2) there are $r$ roots (with multiplicity) of the characteristic equation (15) in $\mathbb{D}$ and $p$ roots (with multiplicity) of the characteristic equation (15) in $\mathcal{U}$.

Remark 4. Under the Cauchy-stability assumption (H1), the inclusion $\Gamma \subset \bar{\mathbb{D}}$ is known. From there, it follows that the unbounded connected component of $\mathbb{C} \setminus \Gamma$ contains the whole set $\mathcal{U}$ so that a weaker form of the lemma is available for considering $z \in \mathcal{U}$ only. If in addition, the considered scheme is also dissipative, meaning that its symbol $\gamma$ satisfies

$$|\gamma(\xi)| \le 1 - \delta|\xi|^{2s}, \quad \xi \in [-\pi, \pi],$$

for some $\delta > 0$ and an integer $s \in \mathbb{N}^*$ independent of $\xi$, then the same separation result is available for $z \in \bar{\mathcal{U}} \setminus \{1\}$. The reason for that property is that in that case one has $\mathbb{S} \cap \Gamma = \{1\}$.

Proof of Lemma 3.

(1) Assume there exists a root $\kappa$ of (15) on the unit circle, then one can find $\theta \in \mathbb{R}$ such that $\kappa = e^{i\theta}$. So we have

$$z = \sum_{j=-r}^{p} a_j \kappa^j = \sum_{j=-r}^{p} a_j e^{ij\theta} = \gamma(\theta).$$

This is a contradiction because $z \in \Gamma$ and by assumption $z \in \mathbb{C} \setminus \Gamma$. This concludes the proof.

(2) We denote $\mathcal{C}$ the unbounded connected component of $\mathbb{C} \setminus \Gamma$. The polynomial (15) has $p+r$ roots (with multiplicity). It is sufficient to count how many roots there are inside the unit disk to deduce the number of roots outside. By continuity of the roots with respect to coefficients and because there is no root on the unit circle for $z \in \mathcal{C}$, we know that there is a constant number of roots inside the unit disk for all $z \in \mathcal{C}$. By Rouché's theorem, one can study the zeros of $f_z(\kappa) = \kappa^r - \frac{1}{z}(a_{-r} + a_{-r+1}\kappa + \cdots + a_p\kappa^{p+r})$ and $g_z(\kappa) = \kappa^r - \frac{1}{z}a_{-r}$ in $\mathbb{D}$ for $z$ sufficiently large to have the result. $\square$

Lemma 3 (Hersh) above is illustrated in Figure 1. The first two lines correspond to the Lemma 3 (Hersh) and the third one describes the possible configuration for $z \in \Gamma \cap \mathbb{S}$, typically not meeting the assumptions. This case will be the object of a subsequent discussion.

For $|z| > 1$, by Lemma 3 (Hersh), the linear subspace of solutions to (14) living in $\ell^2(\mathbb{N})$ is generated by the following $r$ vectors:

$$\begin{pmatrix} 1 \\ \kappa_{\ell} \\ \kappa_{\ell}^2 \\ \kappa_{\ell}^3 \\ \kappa_{\ell}^4 \\ \vdots \end{pmatrix}, \begin{pmatrix} 0 \\ \kappa_{\ell} \\ 2\kappa_{\ell}^2 \\ 3\kappa_{\ell}^3 \\ 4\kappa_{\ell}^4 \\ \vdots \end{pmatrix}, \begin{pmatrix} 0 \\ \kappa_{\ell} \\ 2^2\kappa_{\ell}^2 \\ 3^2\kappa_{\ell}^3 \\ 4^2\kappa_{\ell}^4 \\ \vdots \end{pmatrix}, \dots, \begin{pmatrix} 0 \\ \kappa_{\ell} \\ 2^{\beta_{\ell}-1}\kappa_{\ell}^2 \\ 3^{\beta_{\ell}-1}\kappa_{\ell}^3 \\ 4^{\beta_{\ell}-1}\kappa_{\ell}^4 \\ \vdots \end{pmatrix}, \quad \ell = 1, \dots, M \qquad (16)$$

where $\kappa_1, \dots, \kappa_M$ of multiplicity $\beta_1, \dots, \beta_M$ are the solutions to (15) living in $\mathbb{D}$, with $\beta_1 + \dots + \beta_M = r$ (we omit the $z$-dependence of $\kappa(z)$ for the sake of readability).

**Notation.** We denote $\mathcal{E}^s(z)$ the linear subspace of solutions to (14) living in $\ell^2(\mathbb{N})$ and $K_{i,j}(z) \in \mathcal{M}_{j-i+1,r}(\mathbb{C})$ the matrix where we put in columns the extraction of all the lines between $i$ and $j$ (included) of the $r$ vectors of (16), where $0 \le i \le j$.
or quasi-Toeplitz representation of the scheme. In the same spirit, the computation of the spectral radius of the truncated (i.e. finite dimensional) quasi-Toeplitz matrix may provide significant information for the power boundedness of the method. This is the method used for example by Dakin, Després and Jaouen [11]. This strategy is sometimes called eigenvalue spectrum visualization method, especially by Li, Shu and Zhang [33, 34, 35]. In Section 4, we will compare this latter approach with our own strategy presented hereafter for the O3 scheme in the case of reconstruction boundary conditions. Our strategy is based on the so-called GKS-stability theory introduced by Gustafsson, Kreiss and Sundström [23] which handles the discrete IBVP (1.4)–(1.5)–(1.6) with a zero initial data. The reader can refer to the work by Wu [41] and Coulombel [7] for more recent developments on semi-group estimates in order to deduce the stability of the discrete IBVP (1.4)–(1.5)–(1.6) with non zero initial data from the GKS-stability. The notion of GKS-stability (or also called strong stability) for the boundary problem makes use of the following discrete norms:

$$ \|U_j\|_{\Delta t}^2 = \sum_{n=0}^{+\infty} \Delta t |U_j^n|^2 \quad \text{and} \quad \|U\|_{\Delta x, \Delta t}^2 = \sum_{n=0}^{+\infty} \sum_{j=0}^{+\infty} \Delta t \Delta x |U_j^n|^2. $$

The so-called strong stability, or GKS-stability, is defined by:

**DEFINITION 1.3 (Strong stability).** — The scheme (1.4)–(1.5)–(1.6) is strongly stable if, for ($f_j$) = 0, there exist $C > 0$ and $\alpha_0$, such that for all $\alpha > \alpha_0$, for all boundary data ($g_j^n$), for all $\Delta x > 0$, for all $n \in \mathbb{N}$, the solution satisfies

$$ \sum_{j=-r}^{-1} \|e^{-\alpha n \Delta t} U_j\|_{\Delta t}^2 + \left( \frac{\alpha - \alpha_0}{\alpha \Delta t + 1} \right) \|e^{-\alpha n \Delta t} U\|_{\Delta x, \Delta t}^2 \\ \le C \sum_{j=-r}^{-1} \|e^{-\alpha n \Delta t} g_j\|_{\Delta t}^2. \quad (1.9) $$

We warn the reader that $\|e^{-\alpha n \Delta t} U_j\|_{\Delta t}^2$ is here an abuse of notation to describe $\sum_{n=0}^{+\infty} \Delta t e^{-2\alpha n \Delta t} |U_j^n|^2$ and similarly for $\|e^{-\alpha n \Delta t} U\|_{\Delta x, \Delta t}^2$. This stability definition admits a similar but continuous form for the solutions to continuous hyperbolic PDE's [4]. Namely, it provides some a priori estimates that are useful for a general analysis of such problems.

The following Kreiss theorem [28] expresses a necessary and sufficient condition for the strong stability. We provide hereafter a condensed formulation of this theorem, obtained from [23, Thm. 5.1] combined with [22, Lem. 13.1.4] or with [21, Def. 2.23].
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Illustration of Lemma 3: case $|z| > 1$ (first line), case $|z| = 1$ and $z \notin \Gamma$ (second line) and case $z \in \Gamma$ where Lemma 3 does not hold (third line).</figcaption></figure>

*Remark 5.* For $r = 2$, if the solutions to (15) are $\kappa_1(z) \neq \kappa_2(z)$, then there are exactly two roots with multiplicity 1. The solutions to (14) can be written $\tilde{U}_j(z) = \alpha_1\kappa_1(z)^j + \alpha_2\kappa_2(z)^j$, and we have

$$K_{0,2}(z) = \begin{pmatrix} 1 & 1 \\ \kappa_1(z) & \kappa_2(z) \\ \kappa_1(z)^2 & \kappa_2(z)^2 \end{pmatrix}.$$

*Remark 6.* Still for $r = 2$, if the solution to (15) now is $\kappa(z)$ with multiplicity 2, then the solutions to (14) can be written $\tilde{U}_j(z) = (\alpha_1 + \alpha_2 j)\kappa(z)^j$, and we have

$$K_{0,3}(z) = \begin{pmatrix} 1 & 0 \\ \kappa(z) & \kappa(z) \\ \kappa(z)^2 & 2\kappa(z)^2 \\ \kappa(z)^3 & 3\kappa(z)^3 \end{pmatrix}.$$

We raise awareness of the dependence in $z$ and of the continuity issues because the map $z \mapsto K_{i,j}(z)$ is not continuous whereas the set of roots of (15) is a continuous mapping with respect to $z$.
**THEOREM 1.4** (Kreiss). — The following statements are equivalent:

(i) The scheme (1.4)-(1.5)-(1.6) is strongly stable in the sense of Definition 1.3.

(ii) The Uniform Kreiss–Lopatinskii Condition is satisfied.

The Uniform Kreiss–Lopatinskii Condition corresponds to the absence of zeros for the so-called Kreiss–Lopatinskii determinant $\Delta_{KL}$ that we present here by the informal definition:

$$ \Delta_{KL}(z) = \det(\mathfrak{B}e_1(z), \dots, \mathfrak{B}e_r(z)) \quad (1.10) $$

where $(e_1(z), \dots, e_r(z))$ is an explicit basis of the linear space of the $\ell^2(\mathbb{N})$-stable solutions of the $\mathcal{Z}$-transform of the interior equation (1.4) and $\mathfrak{B}$ is an encoding of the $\mathcal{Z}$-transform of the boundary equation (1.5). For a proper definition of this determinant, the reader can look at Definition 2.9 or the book by Gustafsson, Kreiss and Oliger [22]. Before going on, let us provide some comments to a particular case we already studied.

## 1.2. The case of totally upwind schemes and summary of [5]

The present article is a non trivial extension of our previous work [5] that deals with the restricted case of totally upwind schemes. Totally upwind schemes are schemes of the form (1.2) with $p = 0$ if $a > 0$ or $r = 0$ if $a < 0$. Without loss of generality, we restrict here the discussion to the case $p = 0$ since flipping the indices may turn a case to the other. In this section, we summarize the result of [5] and introduce the novelty of the present work. The first step of the analysis conducted in [5] is based on the introduction of the intrinsic Kreiss–Lopatinskii determinant:

$$ \Delta(z) = \frac{\det(\mathfrak{B}e_1(z), \dots, \mathfrak{B}e_r(z))}{\det(e_1(z), \dots, e_r(z))} \quad (1.11) $$

using the same informal notation as in (1.10). Under appropriate assumptions, an explicit formula for the intrinsic Kreiss–Lopatinskii determinant is obtained:

$$ \forall |z| \ge 1, \quad \Delta(z) = (-1)^{r(m-r)} \det C(z) \left( \frac{a_r}{a_0 - z} \right)^{m-r} \quad (1.12) $$

where $\det C(z)$ is a computable polynomial in $z$ depending only on the coefficients $(a_j)_{j=-r}^0$ and on $\mathfrak{B}$. Thanks to this result, we prove that $\Delta$ is holomorphic on $\bar{U}$. Note that this property may be wrong as long as the standard Kreiss–Lopatinskii determinant is concerned.
Indeed, the root curves $(\kappa_j(z))_j$ can intersect, when a multiple root occurs. For example, for $r = 2$, if there is $(z_n)_{n \in \mathbb{N}} \subset \mathcal{U}$ with $\kappa_1(z_n) \neq \kappa_2(z_n)$ which converge to $z_\infty \in \mathcal{U}$ such that $\kappa_1(z_\infty) = \kappa_2(z_\infty)$ a double root, then we have

$$ \forall j \in \{1, 2\}, \quad \kappa_j(z_n) \xrightarrow[n \to \infty]{} \kappa_j(z_\infty) $$

but

$$ K_{0,3}(z_n) = \begin{pmatrix} 1 & 1 \\ \kappa_1(z_n) & \kappa_2(z_n) \\ \kappa_1^2(z_n) & \kappa_2^2(z_n) \\ \kappa_1^3(z_n) & \kappa_2^3(z_n) \end{pmatrix} \xrightarrow[n \to \infty]{} K_{0,3}(z_\infty) = \begin{pmatrix} 1 & 0 \\ \kappa_1(z_\infty) & \kappa_1(z_\infty) \\ \kappa_1^2(z_\infty) & 2\kappa_1^2(z_\infty) \\ \kappa_1^3(z_\infty) & 3\kappa_1^3(z_\infty) \end{pmatrix}. $$

Consequently, the considered basis (16) of $\mathcal{E}^s(z)$ does not generally define a continuous mapping with respect to $z$.

In spite of the difficulty enlightened above, it turns out that $\mathcal{E}^s(z)$ is a continuous and even holomorphic vector bundle over $\mathcal{U}$ as it is discussed in [6, Thm 4.3]. It is also proved that this vector bundle $\mathcal{E}^s(z)$ can even be continuously extended over $\bar{\mathcal{U}}$, thus considering $z \in \mathbb{S}$ as well (see also [30] for a similar property for the hyperbolic-parabolic PDE case). The main point therein is that for some $z_0 \in \mathbb{S}$, there may exists one (or several) root $\kappa_0(z_0)$ of (15) on $\mathbb{S}$. At such points $z_0$ the Lemma 3 (Hersh) does not hold anymore. This situation is depicted on the third line of Figure 1. For $z$ on $\mathbb{S}$, the space $\mathcal{E}^s(z)$ still is of dimension $r$ and we extend the notation $K_{i,j}(z)$. We can summarize the above discussion in the following theorem.

**Theorem 7** ([6]). Under assumption (H1), the space $\mathcal{E}^s(z)$ is a holomorphic vector bundle over $\mathcal{U}$ and can be extended in a unique way to a continuous vector bundle over $\bar{\mathcal{U}}$.

**Remark 8.** For the extension, the first difficulty is to select the roots of (15) coming from the inside, indeed, if there is a root on $\mathbb{S}$, it can be coming from the inside of $\mathbb{D}$, the outside or both (in case of multiplicity). In Section 4.2, we will explain the numerical strategy to select the good ones. The second difficulty is to prove the continuity of $\mathcal{E}^s(z)$ after the extension, it follows from the existence of a K-symmetrizer and is obtained e.g. in [6, Thm 4.3]. As previously observed, $K_{i,j}(z)$ is generally not continuous with respect to $z$.

## 2.2. Intrinsic Kreiss-Lopatinskii determinant.
In this section, we define properly formulas (10) and (11). Let us consider the $\mathcal{Z}$-transformed version of the boundary condition (7), that is

$$ z \begin{pmatrix} \tilde{U}_0(z) \\ \vdots \\ \tilde{U}_{r-1}(z) \end{pmatrix} - \mathcal{B} \begin{pmatrix} \tilde{U}_0(z) \\ \vdots \\ \tilde{U}_{m-1}(z) \end{pmatrix} = \begin{pmatrix} a_{-r} & \cdots & a_{-1} \\ \vdots & \ddots & \vdots \\ 0 & a_{-r} \end{pmatrix} \begin{pmatrix} \tilde{g}_{-r}(z) \\ \vdots \\ \tilde{g}_{-1}(z) \end{pmatrix}. \quad (17) $$

Injecting the solution $(\tilde{U}_j(z))_{j \in \mathbb{N}} \in \mathcal{E}^s(z)$ to (14) into (17), we obtain a system of $r$ equations with $r$ scalar unknowns: they are the coefficients of $(\tilde{U}_j(z))_{j \in \mathbb{N}}$ written in the basis (16) of $\mathcal{E}^s(z)$.

**Remark 9.** For $r = 2$ and a given value of $z$ (we skip for convenience the dependence in $z$ hereafter), if $\kappa_1 \neq \kappa_2$ so that the solution to (17) has the form $\alpha_1 \kappa_1^j + \alpha_2 \kappa_2^j$, then that solution is constrained by the system (17). The matricial form of that system reads

$$ z \begin{pmatrix} 1 & 1 \\ \kappa_1 & \kappa_2 \\ \kappa_1^2 & \kappa_2^2 \\ \vdots & \vdots \\ \kappa_1^{m-1} & \kappa_2^{m-1} \end{pmatrix} - \mathcal{B} \begin{pmatrix} 1 \\ \kappa_1 \\ \kappa_1^2 \\ \vdots \\ \kappa_1^{m-1} \end{pmatrix} = \begin{pmatrix} \alpha_1 \\ \alpha_2 \end{pmatrix} = \begin{pmatrix} a_{-1} \tilde{g}_{-1} + a_{-2} \tilde{g}_{-2} \\ a_{-2} \tilde{g}_{-1} \end{pmatrix}. $$

The injectivity, whence invertibility, of the boundary condition is thus directly related to the property $\det(zK_{0,1} - BK_{0,m-1}(z)) \neq 0$, where $zK_{0,1} - BK_{0,m-1}(z) \in \mathcal{M}_{2,2}(\mathbb{C})$.
Applying the residue theorem to $\Delta$, we developed a numerical strategy to count the number of zeros of the Kreiss–Lopatinskii determinant in $U$. By Theorem 1.4 (Kreiss), we conclude that if $\text{Ind}_{\Delta(S)}(0) < r$ then the scheme is not stable where $\text{Ind}_{\Delta(S)}(0)$ is the notation for the winding number of 0 with respect to the Kreiss–Lopatinskii curve $\Delta(S)$. This result allows us to establish an efficient and practical method (see Method 19 of [5] or Method 2.14 of the present paper) to study the stability of a scheme with boundary. It provides sharp results for the solution $(U_j^n)_{n \in \mathbb{N}} \in \ell^2(\mathbb{N})$ to the problem set on the half line $\mathbb{N}$. In particular, contrary to numerical investigations of stability which are based on the computation of the spectral radius, no arbitrary truncation of (quasi-)Toeplitz matrices is needed. In return, a problem set on a bounded space domain needs, for a whole convergence study, superposition techniques for truncated data, as used in [6] and [8]. This feature mainly restricts the study to explicit scheme.

In the present article, we extend to $p \ge 1$ the previous results of [5]. The main improvement therein is now to cover schemes with arbitrary high order. As a consequence of the Iserles–Strang theorem [27], the order of accuracy of totally upwind schemes is restricted to 2. Using appropriate techniques such as the Lax–Wendroff procedure, also known as the Cauchy–Kowalewski procedure, high order both in space and time can be reached even with one-step methods, see for instance [12, 36]. We prove that the stability of the scheme is, once again, fully characterized by the winding number of 0, an easily computable quantity. Even if an explicit formula such as (1.12) is now generally unavailable, the holomorphic property of the intrinsic Kreiss–Lopatinskii determinant remains true and sufficies to build the same method for studying the stability. This approach, using the winding number, is robust since instead of finding zeros of an algebraic curve, it only requires the computation of a winding number, which is an integer, to count the number of zeros. One can mention the work of Thuné [39] who develops a numerical method to check the GKS-stability. He looks for the precise location of the zeros of the Kreiss–Lopatinskii determinant approximating the roots of some parameterized characteristic polynomial equations which is significantly different with our work. One can also cite the work of Tadmor and Goldberg [19, 20] which gives a sufficient condition to have the stability of a scheme with boundaries. One can see this condition as a weaker version of the Kreiss–Lopatinskii determinant but it has the advantage of decorrelating the study of the interior equation and the boundary equations. In our paper, the framework is less restrictive since the scheme is not necessary dissipative and the boundary conditions is more general.
**Definition 10** (Kreiss-Lopatinskii determinant). The Kreiss-Lopatinskii determinant is the complex-valued function defined for $|z| \ge 1$ by:

$$ \Delta_{\text{KL}}(z) = \det(zK_{0,r-1}(z) - BK_{0,m-1}(z)). $$

Despite the fact that the space $\mathcal{E}^s(z)$ is a holomorphic vector bundle over $U$ and continuous over $\bar{U}$ (Theorem 7), this determinant $\Delta_{\text{KL}}$ is not holomorphic on $U$. To retrieve those properties, we define the intrinsic Kreiss-Lopatinskii determinant $\Delta$ that we can motivate by the following informal discussion. The above Kreiss-Lopatinskii determinant is actually not well defined until we order in some way the roots $(\kappa_j(z))_{j=1,...,r}$ of (15). There are two points to emphasize. The first one is related to crossing roots and already discussed after Remark 6. The second one is that, outside crossing cases, being given any choice for the ordering of the roots (and thus of the vectors of the basis (16) for the vector bundle), there is in general no chance to obtain a holomorphicity property for the components of the matrix $K_{0,m-1}(z)$ over $U$. For example, even the roots of $X^2 - z$ are not holomorphic w.r.t $z \in U$ because of the logarithm determination. On the other side, any symmetric functions of the roots $(\kappa_j(z))_{j=1,...,r}$ however are holomorphic because they can be obtained directly in terms of the coefficients of the polynomial (15). So, except for crossing roots, the same holds for the quantity $\Delta_{\text{KL}}(z)$ since the matrices $B$ and $\tilde{B}$ are constants and the determinant itself is a symmetric function.

A very natural way to reach the holomorphic property and go beyond the last difficulties consists in dividing $\Delta_{\text{KL}}$ by the quantity $\det K_{0,r-1}(z)$. Hence, the same permutation or combination of the vectors of the basis (16) is involved in both computations. This intrinsic Kreiss-Lopatinskii determinant has already been introduced and studied in a particular case in [4].

**Definition 11** (Intrinsic Kreiss-Lopatinskii determinant). The intrinsic Kreiss-Lopatinskii determinant is the complex-valued function defined for $|z| \ge 1$ by:

$$ \Delta(z) = \frac{\Delta_{\text{KL}}(z)}{\det K_{0,r-1}(z)}. \qquad (18) $$

Let us note that the intrinsic Kreiss-Lopatinskii determinant can be rewritten

$$ \Delta(z) = \frac{\det(zK_{0,r-1}(z) - BK_{0,m-1}(z))}{\det K_{0,r-1}(z)} = z^r \det \left( I_r - \frac{BK_{0,m-1}(z)K_{0,r-1}(z)^{-1}}{z} \right). \quad (19) $$

To conclude with these definitions, let us state a little more about the Uniform Kreiss-Lopatinskii Condition. With the above notations and additionally to the invertibility of $zK_{0,r-1}(z) - BK_{0,m-1}(z)$, it corresponds to the existence of a constant $C > 0$ such that for any $z \in \bar{U}$, any $\tilde{U} \in \mathcal{E}^s(z)$ solution to (17) satisfies the uniform estimate $\|\tilde{U}\| \le C\|\tilde{g}\|$. From the Parseval identity for the $Z$-transform, this inequality gives directly the first necessary half-part of the strong stability estimate (9). We refer the reader to [18] for a more detailed presentation.

**2.3. Main results.** Theorem 12 is our main theoretical result. It states that the intrinsic Kreiss-Lopatinskii determinant has the same regularity properties as $\mathcal{E}^s(z)$, see Theorem 7.

**Theorem 12** (Smoothness of the intrinsic Kreiss-Lopatinskii determinant). Assume (H1). The intrinsic Kreiss-Lopatinskii determinant $\Delta$ is holomorphic on $U$ and continuous on $\bar{U}$.

By equation (18), the function $\Delta$ shares the same zeros with the Kreiss-Lopatinskii determinant $\Delta_{\text{KL}}$, so that it can be used as an alternative in the Uniform Kreiss-Lopatinskii Condition, see Theorem 2 (Kreiss). Another property, important for the forthcoming applications, lies in the next Corollary 14 and involves the following important geometrical object:

**Definition 13.** The Kreiss-Lopatinskii curve $\Delta(\mathcal{S})$ is the closed complex parameterized curve

$$ \Delta(\mathcal{S}) = \{\theta \in [0, 2\pi] \mapsto \Delta(e^{i\theta})\}. $$
## 1.3. Outline of the paper

After constructing the intrinsic Kreiss–Lopatinskii determinant in Section 2 and Section 3, we see that, in such a general case, the lack of an explicit formula for $\Delta(z)$ does not preclude holomorphic properties (see Theorem 2.11). From there, we obtain the following stability criterion: if $\mathrm{Ind}_{\Delta(S)}(0) < r$ then the scheme is not stable (see Corollary 2.13). We prove these results by the use of Hermite interpolation and residue theorem. To compute the Kreiss–Lopatinskii determinant numerically, in Section 4, we use an easy-to-use formulation of it which is, in some sense, close to the explicit formulation (1.12). Moreover, Section 4 gathers the numerical procedure to draw the Kreiss–Lopatinskii curve, several examples and numerical experiments for illustrating the efficiency of the proposed strategy.

## 2. Kreiss–Lopatinskii determinants

In this section, we introduce the Kreiss–Lopatinskii determinant, a usual tool to check the Uniform Kreiss–Lopatinskii Condition. Then we define the intrinsic Kreiss–Lopatinskii determinant, namely a reshaping of the previous one, which is more convenient in practice since it enjoys improved properties compared to the regular Kreiss–Lopatinskii determinant: holomorphicity, continuity, independence on the basis...

### 2.1. Stable subspace $\mathcal{E}^s(z)$ and matrix representation

First, we study the solutions to the interior equation:

$$U_j^{n+1} = \sum_{k=-r}^{p} a_k U_{k+j}^n, \quad j \in \mathbb{N}, n \in \mathbb{N}. \qquad (2.1)$$

To study this equation, the $\mathcal{Z}$-transform (see [17, Lesson 40]) is applied. This transformation is defined for $(x_n)_{n \in \mathbb{N}} \in \ell^2(\mathbb{N})$ such that $x_0 = 0$ and $z \in \mathcal{U}$ by $\tilde{x}(z) = \sum_{n \ge 0} z^{-n} x_n$. The previous equation then reads

$$z \tilde{U}_j(z) = \sum_{k=-r}^{p} a_k \tilde{U}_{j+k}(z), \quad j \in \mathbb{N}, z \in \mathcal{U}. \qquad (2.2)$$

To solve the linear recurrence equation (2.2), let us introduce the following characteristic equation where $z$ plays the role of a parameter and $\kappa$ is the indeterminate:

$$z\kappa^r = \sum_{k=-r}^{p} a_k \kappa^{r+k}. \qquad (2.3)$$
Using the residue theorem<fnref n="1" /> thanks to Theorem 12, we obtain the following result.

**Corollary 14** (Number of zeros of the intrinsic Kreiss-Lopatinskii determinant). Assume (H1). If $0 \notin \Delta(\mathbb{S})$ then the equation $\Delta(z) = 0$ has exactly $r - \text{Ind}_{\Delta(\mathbb{S})}(0)$ zeros in $\mathcal{U}$.

Here above and in all the paper, $\text{Ind}_{\Delta(\mathbb{S})}(0)$ denotes the winding number of the origin with respect to the closed oriented curve $\Delta(\mathbb{S})$ (see [24] for a definition of the winding number). This corollary helps us to establish an efficient and practical method to study the stability of a given IBVP through Theorem 2 (Kreiss). In particular, the low computational cost of the following procedure is very appealing for the study of parameterised IBVP's, see Section 4.

**Method 15** (Uniform Kreiss-Lopatinskii Condition check). There are two different cases:

* if $0 \in \Delta(\mathbb{S})$, then there exists $z_0 \in \mathbb{S}$ such that $\Delta(z_0) = 0$.

* if $0 \notin \Delta(\mathbb{S})$, $\Delta$ does not vanish on $\mathbb{S}$ and it has $r - \text{Ind}_{\Delta(\mathbb{S})}(0)$ zeros in $\mathcal{U}$ by Corollary 14. It follows that if $\text{Ind}_{\Delta(\mathbb{S})}(0) = r$ then the scheme is stable. Otherwise the scheme is unstable.

In summary, by Theorem 2 (Kreiss) and since Uniform Kreiss-Lopatinskii Condition is fulfilled if and only if the Kreiss-Lopatinskii determinant has no zero in $\mathcal{U}$, Method 15 can be used to conclude that the scheme is stable or not. Some illustrations for the O3 scheme and the fifth-order Lax-Wendroff scheme follow in Section 4.

### 3. PROOF OF THEOREM 12 AND COROLLARY 14

**3.1. Constant-recursive sequence of order r.** For each $z \in \mathcal{U}$, we denote $P_z$ the polynomial linked to the characteristic equation (15), i.e.

$$P_z(\kappa) = a_p\kappa^{r+p} + \dots + a_1\kappa^{r+1} + (a_0 - z)\kappa^r + a_{-1}\kappa^{r-1} + \dots + a_{-r+1}\kappa + a_{-r}. \quad (20)$$

By Lemma 3 (Hersh), the polynomial $P_z(\kappa)$ can be factorized into two polynomials: one with the $r$ roots in $\mathbb{D}$, denoted $R_z(\kappa)$ and one with the $p$ roots in $\mathcal{U}$, denoted $Q_z(\kappa)$. We know that the coefficients of $P_z$ are holomorphic in $z$. We already said that the basis (16) is not holomorphic because the roots $\kappa$ are not. In the next result, we prove that the symmetric functions of the $r$ roots $\kappa$ living in $\mathbb{D}$ are indeed holomorphic in $\mathcal{U}$, in other words, the coefficients of $R_z$ are holomorphic in $\mathcal{U}$.

**Lemma 16.** For all $z \in \mathcal{U}$, the polynomial $R_z(X) = \prod_{j=1}^r (X - \kappa_j(z))$ has holomorphic coefficients in $\mathcal{U}$, where $(\kappa_j(z))_{j=1}^r$ are the $r$ roots (with multiplicity) in $\mathbb{D}$ of (15).

*Proof.* We use the Dunford-Taylor formula with $C_z$ the companion matrix of the polynomial (15):

$$\Pi(z) = \frac{1}{2\pi} \int_{\mathbb{S}} (\zeta I_{r+p} - C_z)^{-1} d\zeta$$

It is the projection along $\mathbb{E}_s(z) = \ker \prod_{j=1}^r (C_z - \kappa_j(z))$ onto $\mathbb{E}_u(z) = \ker \prod_{j=r+1}^{r+p} (C_z - \kappa_j(z))$ where $(\kappa_j(z))_{j=r+1}^{r+p}$ are the roots of (15) in $\mathcal{U}$, because $(\kappa_j(z))_{j=1}^r$ are surrounded by $\mathbb{S}$ and $(\kappa_j(z))_{j=r+1}^{r+p}$ are not. The projector $\Pi(z)$ is holomorphic on $\mathcal{U}$ since it is a holomorphic parameter integral. We have $C_z \circ \Pi(z)|_{\mathbb{E}_u(z)} = 0$ and $C_z \circ \Pi(z)|_{\mathbb{E}_s(z)} = C_z$, then the characteristic polynomial of $C_z \circ \Pi(z)$ is $X^p R_z(X)$ because $\mathbb{C}^{r+p} = \mathbb{E}_s(z) \oplus \mathbb{E}_u(z)$. The function $z \mapsto C_z \circ \Pi(z)$ is holomorphic on $\mathcal{U}$, then the coefficients of its characteristic polynomial are too. This concludes the proof. $\square$

<footnote n="1">¹All the complex analysis results can be found in [24].</footnote>
This equation is nothing but the discrete dispersion relation of the finite difference scheme (2.1), with frequency parameter $κ$ in space and $z$ in time. It is formally obtained by looking for solutions to the interior equation (2.1) having the form $U_j^n = z^n \kappa^j$.

In the spirit of a common result by Hersh [26], the following lemma provides a splitting result for the roots with respect to the unit circle.

**LEMMA 2.1 (Hersh).** — Assume 1.2. For $z$ in the unbounded connected component of $C \setminus \Gamma$,

(1) there is no root of the characteristic equation (2.3) on $\mathbb{S}$,

(2) there are $r$ roots (with multiplicity) of the characteristic equation (2.3) in $\mathbb{D}$ and $p$ roots (with multiplicity) of the characteristic equation (2.3) in $U$.

**Remark 2.2.** — Under the Cauchy-stability assumption 1.2, the inclusion $\Gamma \subset \bar{\mathbb{D}}$ is known. From there, it follows that the unbounded connected component of $C \setminus \Gamma$ contains the whole set $U$ so that a weaker form of the lemma is available for considering $z \in U$ only. If in addition, the considered scheme is also dissipative, meaning that its symbol $\gamma$ satisfies

$$|\gamma(\xi)| \le 1 - \delta|\xi|^{2s}, \quad \xi \in [-\pi, \pi],$$

for some $\delta > 0$ and an integer $s \in \mathbb{N}^*$ independent of $\xi$, then the same separation result is available for $z \in \bar{U} \setminus \{1\}$. The reason for that property is that in that case one has $\mathbb{S} \cap \Gamma = \{1\}$.

*Proof of Lemma 2.1.* —

(1) Assume there exists a root $\kappa$ of (2.3) on the unit circle, then one can find $\theta \in \mathbb{R}$ such that $\kappa = e^{i\theta}$. So we have

$$z = \sum_{j=-r}^{p} a_j \kappa^j = \sum_{j=-r}^{p} a_j e^{ij\theta} = \gamma(\theta).$$

This is a contradiction because $z \in \Gamma$ and by assumption $z \in C \setminus \Gamma$. It concludes the proof.

(2) We denote $\mathcal{C}$ the unbounded connected component of $C \setminus \Gamma$. The polynomial (2.3) has $p+r$ roots (with multiplicity). It is sufficient to count how many roots there are inside the unit disk to deduce the number of roots outside. By continuity of the roots with respect to coefficients and because there is no root on the unit circle for $z \in \mathcal{C}$, we know that there is a constant number of roots inside the unit disk for all $z \in \mathcal{C}$. By Rouché's theorem, one can study the zeros of $f_z(\kappa) = \kappa^r - \frac{1}{z}(a_{-r} + a_{-r+1}\kappa + \cdots + a_p\kappa^{p+r})$ and $g_z(\kappa) = \kappa^r - \frac{1}{z}a_{-r}$ in $\mathbb{D}$ for $z$ sufficiently large to have the result. $\square$
**3.2. Hermite interpolation.** To prove the holomorphic properties of $\Delta$, by (19), it is sufficient to study the function $z \mapsto K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$. To simplify this study, for $z \in \bar{\mathcal{U}}$, we introduce on the linear map

$$ \varphi_z : Q \in \mathbb{C}_{m-1}[X] \mapsto \varphi_z(Q) \in \mathbb{C}_{r-1}[X] \quad (21) $$

where $\varphi_z(Q)$ is the Hermite interpolation polynomial of degree less than $r-1$ defined by the value $(Q(\kappa_1(z)), \dots, Q(\beta_{1-1})(\kappa_1(z)), Q(\kappa_2(z)), \dots, Q(\beta_{2-1})(\kappa_2(z)), \dots, Q(\beta_{M-1})(\kappa_M(z)))$, where the $\kappa$'s are the same as in (16). The link between the matrix $K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$ and $\varphi_z$ is given by:

**Lemma 17.** For all $z \in \bar{\mathcal{U}}$, the transpose of the matrix $K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$ is the representation in the canonical basis of $\varphi_z$ defined in (21).

*Proof.* The Hermite interpolation make appear the following matrix

$$ H_{0,j}(z) = \underbrace{\begin{pmatrix} 1 & 0 & 0 & \cdots & 1 & 0 & \cdots & 1 & 0 & \cdots \\ \kappa_1 & 1 & 0 & \cdots & \kappa_2 & 1 & \cdots & \kappa_M & 1 & \cdots \\ \kappa_1^2 & 2\kappa_1 & 2 & \cdots & \kappa_2^2 & 2\kappa_2 & \cdots & \kappa_M^2 & 2\kappa_M & \cdots \\ \kappa_1^3 & 3\kappa_1^2 & 6\kappa_1 & \cdots & \kappa_2^3 & 3\kappa_2^2 & \cdots & \kappa_M^3 & 3\kappa_M^2 & \cdots \\ \kappa_1^4 & 4\kappa_1^3 & 12\kappa_1^2 & \cdots & \kappa_2^4 & 4\kappa_2^3 & \cdots & \kappa_M^4 & 4\kappa_M^3 & \cdots \\ \vdots & \vdots & \vdots & \cdots & \vdots & \vdots & \cdots & \vdots & \vdots & \cdots \\ \kappa_1^j & j\kappa_1^{j-1} & j(j-1)\kappa_1^{j-2} & \cdots & \kappa_2^j & j\kappa_2^{j-1} & \cdots & \kappa_M^j & j\kappa_M^{j-1} & \cdots \end{pmatrix}}_{\beta_1 \text{ columns linked to } \kappa_1} \underbrace{\begin{pmatrix} \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \end{pmatrix}}_{\beta_2 \text{ columns linked to } \kappa_2} \underbrace{\begin{pmatrix} \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \end{pmatrix}}_{\beta_M \text{ columns linked to } \kappa_M} $$

The representation of $z \mapsto \varphi_z$ in the canonical basis is $(H_{0,m-1}(z)H_{0,r-1}^{-1}(z))^T$. Besides, there exists an invertible matrix $M(z) \in \mathcal{M}_r(\mathbb{C})$ such that $K_{0,j}(z) = H_{0,j}(z)M(z)$. Therefore, we have $K_{0,m-1}(z)K_{0,r-1}^{-1}(z) = H_{0,m-1}(z)M(z)M(z)^{-1}H_{0,r-1}^{-1}(z) = H_{0,m-1}(z)H_{0,r-1}^{-1}(z)$. The result follows. $\square$

**Proposition 18.** The function $z \mapsto \varphi_z$ is holomorphic on $\mathcal{U}$.

*Proof.* For all $k \in [0 : m-1]$, we want every coefficient of the polynomial $\varphi_z(X^k)$ to be holomorphic on $\mathcal{U}$. Writing $\varphi_z(X^k)(x) = \sum_{j=0}^{r-1} \alpha_{j,k}(z)x^j$, we know that

$$ \forall j \in [0 : r-1], \quad j!\alpha_{j,k}(z) = \partial_x^j \varphi_z(X^k)(x)|_{x=0}. \quad (22) $$

By the error of Hermite interpolation (see [21]), we have

$$ \varphi_z(X^k)(x) - x^k = \frac{1}{2i\pi} \int_{\mathcal{S}} \frac{\zeta^k R_z(x)}{(x-\zeta)R_z(\zeta)} d\zeta \quad (23) $$

where $R_z(X)$ is defined in Lemma 16. Differentiating equation (23) (with the Leibniz product rule), one obtains

$$ j!\alpha_{j,k}(z) = k!\delta_k^j + \sum_{s=0}^{j} \binom{j}{s} R_z^{(j-s)}(0) \frac{1}{2i\pi} \int_{\mathcal{S}} \frac{-s! \zeta^{k-s-1}}{R_z(\zeta)} d\zeta. \quad (24) $$

By Lemma 16 and the holomorphicity of parameter-dependent integrals, the function $z \mapsto \alpha_{j,k}(z)$ is holomorphic on $\mathcal{U}$ for all $j \in [0 : r-1]$ and $k \in [0 : m-1]$. The proof is now complete. $\square$

**Proposition 19.** The function $z \mapsto \varphi_z$ is continuous on $\bar{\mathcal{U}}$.

*Proof.* Because Lemma 3 (Hersh) does not hold anymore for $z \in \mathcal{S}$, the roots $\kappa(z)$ of characteristic equation (15) can be on the unit circle $\mathcal{S}$. To prove the continuity of $z \mapsto \alpha_{j,k}(z)$, we use equation (24) but replacing $\mathcal{S}$ by $\mathcal{S}_\varepsilon \stackrel{\text{def}}{=} \{z \in \mathbb{C}, |z| = 1+\varepsilon\}$ for $\varepsilon > 0$. Using the continuity of parameter-dependent integrals and the continuity of the roots $\kappa(z)$ of characteristic equation (15), we obtain the continuity
<figure><img src="image_5.png" /><figcaption>Figure 2.1. Illustration of Lemma 2.1: case $|z| > 1$ (first line), case $|z| = 1$ and $z \notin \Gamma$ (second line) and case $z \in \Gamma$ where Lemma 2.1 does not hold (third line).</figcaption></figure>

Lemma 2.1 (Hersh) above is illustrated in Figure 2.1. The first two lines correspond to the Lemma 2.1 (Hersh) and the third one describes the possible configuration for $z \in \Gamma \cap \mathcal{S}$, typically not meeting the assumptions. This case will be the object of a subsequent discussion.
of the coefficients of $R_z$ and thus the function $z \mapsto \alpha_{j,k}(z)$ is continuous on $\bar{\mathcal{U}}$ for all $j \in [0 : r - 1]$ and $k \in [0 : m - 1]$. The proof is now complete. $\square$

**Proposition 20.** The function $z \mapsto \varphi_z$ is bounded on $\bar{\mathcal{U}}$.

*Proof.* Equation (24) can give a bound of every components of $K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$. Indeed, using Rouché's theorem, as in the proof of Lemma 3 (Hersh), we can see that for $|z| > R$ for a certain $R$, all the roots $\kappa(z)$ of the characteristic equation (15) satisfy $|\kappa(z)| < \frac{1}{2}$. Then, for $|z| > R$, one can have

$$ \left| \frac{1}{2i\pi} \int_{\mathcal{S}} \frac{-s! \zeta^{k-s-1}}{R_z(\zeta)} d\zeta \right| = \left| \frac{s!}{2i\pi} \int_0^{2\pi} \frac{e^{i\theta(k-s-1)}}{\prod_{j=1}^r (e^{i\theta} - \kappa_j(z))} d\theta \right| \le \frac{s!}{2\pi} \int_0^{2\pi} \frac{1}{|1 - \frac{1}{2}|^r} d\theta \le s! 2^r. $$

By Gauss-Lucas theorem, the roots of all the derivatives of $R_z$ are in $\mathbb{D}$ for all $|z| > 1$, it follows that $R_z^{(j-s)}(0)$ is bounded independently of $z$. Then for $|z| > R$, the quantity $K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$ is bounded. Moreover, by Proposition 19, the quantity $K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$ is bounded on the compact set $\{1 \le |z| \le R\}$. The proof is now complete. $\square$

### 3.3. Conclusion.

*Proof of Theorem 12.* By Lemma 17, the continuity and holomorphicity properties of $z \mapsto \varphi_z$ provided in Propositions 18 and 19 are shared by the function $z \mapsto K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$. The expression of the intrinsic determinant (19) concludes the proof. $\square$

The next proof is close to the proof of the Corollary 15 of [4]. We reproduce it here for completeness.

*Proof of Corollary 14.* Let us define the following function

$$ \tilde{\Delta} : z \in \mathbb{D}^* \mapsto \Delta(1/z) \in \mathbb{C}. $$

By Theorem 12, the function $\tilde{\Delta}$ is meromorphic on $\mathbb{D}$ with one only pole in 0 and is continuous on $\bar{\mathbb{D}} \setminus \{0\}$. By Proposition 20, the function $z \mapsto K_{0,m-1}(1/z)K_{0,r-1}(1/z)^{-1}$ is bounded on $\bar{\mathbb{D}}$, it follows that 0 is a pole of order $r$ of the function $\tilde{\Delta}$. The residue theorem applied on $\tilde{\Delta}$ with the path $\mathcal{S}$ gives the following equality:

$$ \mathrm{Ind}_{\tilde{\Delta}(\mathcal{S})}(0) = \#\mathrm{zeros}_{\tilde{\Delta}}(\mathbb{D}) - \#\mathrm{poles}_{\tilde{\Delta}}(\mathbb{D}). $$

It follows that

$$ \#\mathrm{zeros}_{\Delta}(\mathcal{U}) = r - \mathrm{Ind}_{\Delta(\mathcal{S})}(0). $$

$\square$

## 4. NUMERICAL RESULTS

### 4.1. New formulation of $\Delta$.
In [4], an explicit formula of the Kreiss-Lopatinskii determinant is given. Unfortunately to reduce the boundary matrix $\mathcal{B}$, the characteristic equation of degree $r + p$ is used to find a final matrix of size $r \times (r + p)$ and, since in the present analysis $p \ne 0$, the matrix is not square. To skirt that problem we will use the polynomial $R_z$ defined in Lemma 16 instead of using the complete characteristic equation (20). It reads also

$$ R_z(X) = \prod_{j=1}^{r} (X - \kappa_j(z)) = X^r + \sigma_{r-1}(z)X^{r-1} + \cdots + \sigma_1(z)X + \sigma_0(z) $$

where $(\sigma_j(z))_j$ are the symmetric functions of $(\kappa_j(z))_j$.

Because $\tilde{U}_j(z)$ is in $\mathcal{E}^s(z)$ and can be expressed in the basis (16), we have, for all $j \in \mathbb{N}$,

$$ \tilde{U}_{j+r}(z) + \sigma_{r-1}(z)\tilde{U}_{j+r-1}(z) + \cdots + \sigma_1(z)\tilde{U}_{j+1}(z) + \sigma_0(z)\tilde{U}_j(z) = 0. \quad (25) $$
For $|z| > 1$, by Lemma 2.1 (Hersh), the linear subspace of solutions to (2.2) living in $ℓ^2(\mathbb{N})$ is generated by the following $r$ vectors:

$$ \begin{pmatrix} 1 \\ \kappa_{\ell} \\ \kappa_{\ell}^2 \\ \kappa_{\ell}^3 \\ \kappa_{\ell}^4 \\ \vdots \end{pmatrix}, \begin{pmatrix} 0 \\ \kappa_{\ell} \\ 2\kappa_{\ell}^2 \\ 3\kappa_{\ell}^3 \\ 4\kappa_{\ell}^4 \\ \vdots \end{pmatrix}, \begin{pmatrix} 0 \\ \kappa_{\ell} \\ 2^2\kappa_{\ell}^2 \\ 3^2\kappa_{\ell}^3 \\ 4^2\kappa_{\ell}^4 \\ \vdots \end{pmatrix}, \dots, \begin{pmatrix} 0 \\ \kappa_{\ell} \\ 2^{\beta_{\ell}-1}\kappa_{\ell}^2 \\ 3^{\beta_{\ell}-1}\kappa_{\ell}^3 \\ 4^{\beta_{\ell}-1}\kappa_{\ell}^4 \\ \vdots \end{pmatrix}, \quad \ell = 1, \dots, M \quad (2.4) $$

where $\kappa_1, \dots, \kappa_M$ of multiplicity $\beta_1, \dots, \beta_M$ are the solutions to (2.3) living in $\mathbb{D}$, with $\beta_1 + \dots + \beta_M = r$ (we omit the $z$-dependence of $\kappa(z)$ for the sake of readability).

**Notation 2.3.** — We denote $\mathcal{E}^s(z)$ the linear subspace of solutions to (2.2) living in $ℓ^2(\mathbb{N})$ and $K_{i,j}(z) \in \mathcal{M}_{j-i+1,r}(\mathbb{C})$ the matrix where we put in columns the extraction of all the lines between $i$ and $j$ (included) of the $r$ vectors of (2.4), where $0 \le i \le j$.

**Remark 2.4.** — For $r = 2$, if the solutions to (2.3) are $\kappa_1(z) \neq \kappa_2(z)$, then there are exactly two roots with multiplicity 1. The solutions to (2.2) can be written $\tilde{U}_j(z) = \alpha_1\kappa_1(z)^j + \alpha_2\kappa_2(z)^j$, and we have

$$ K_{0,2}(z) = \begin{pmatrix} 1 & 1 \\ \kappa_1(z) & \kappa_2(z) \\ \kappa_1(z)^2 & \kappa_2(z)^2 \end{pmatrix}. $$

**Remark 2.5.** — Still for $r = 2$, if the solution to (2.3) now is $\kappa(z)$ with multiplicity 2, then the solutions to (2.2) can be written $\tilde{U}_j(z) = (\alpha_1 + \alpha_2 j)\kappa(z)^j$, and we have

$$ K_{0,3}(z) = \begin{pmatrix} 1 & 0 \\ \kappa(z) & \kappa(z) \\ \kappa(z)^2 & 2\kappa(z)^2 \\ \kappa(z)^3 & 3\kappa(z)^3 \end{pmatrix}. $$

We raise awareness of the dependence in $z$ and of the continuity issues because the map $z \mapsto K_{i,j}(z)$ is not continuous whereas the set of roots of (2.3) is a continuous mapping with respect to $z$. Indeed, the root curves $(\kappa_j(z))_j$ can intersect, when a multiple root occurs. For example, for $r = 2$, if there is $(z_n)_{n \in \mathbb{N}} \subset \mathcal{U}$ with $\kappa_1(z_n) \neq \kappa_2(z_n)$ which converge to $z_\infty \in \mathcal{U}$ such that $\kappa_1(z_\infty) = \kappa_2(z_\infty)$ a double root, then we have

$$ \forall j \in \{1, 2\}, \quad \kappa_j(z_n) \xrightarrow[n \to \infty]{} \kappa_j(z_\infty) $$
**Notation.** We note, for all $j \in \mathbb{N}$, $\tilde{\mathfrak{U}}_j(z)$ the vector $(\tilde{U}_0(z) \cdots \tilde{U}_j(z))^T$ of size $j+1$.

**Proposition 21.** Let $\mathcal{B} \in \mathcal{M}_{r,m}(\mathbb{C})$. There exists a function $\tilde{\mathcal{B}}: \mathbb{C}^r \to \mathcal{M}_{r,r}(\mathbb{C})$ constructible such that, for all $z \in \bar{\mathcal{U}}$, we have

$$ \mathcal{B}\tilde{\mathfrak{U}}_{m-1}(z) = \tilde{\mathcal{B}}(\sigma_0(z), \dots, \sigma_{r-1}(z))\tilde{\mathfrak{U}}_{r-1}(z) \qquad (26) $$

where $(\tilde{U}_j(z))_j$ satisfies (25) for all $j \in \mathbb{N}$.

By "constructible function", we mean here that we establish a computable algorithm to get the matrix $\tilde{\mathcal{B}}(\sigma_0(z), \dots, \sigma_{r-1}(z))$. This algorithm, based on a Gaussian elimination, is fully described in the following proof.

*Proof.* For $z \in \bar{\mathcal{U}}$ and $\zeta_0 = \sigma_0(z), \dots, \sigma_{r-1} = \sigma_{r-1}(z)$. By a descending induction on $j$ between $m-1$ to $r-1$, we construct a matrix $\mathcal{B}_j(\zeta_0, \dots, \zeta_{r-1}) \in \mathcal{M}_{r,j+1}(\mathbb{C})$ such that

$$ \mathcal{B}\tilde{\mathfrak{U}}_{m-1}(z) = \mathcal{B}_j(\zeta_0, \dots, \zeta_{r-1})\tilde{\mathfrak{U}}_j(z). $$

*Initialization:* if $j = m-1$ then one can take $\mathcal{B}$ for the matrix $\mathcal{B}_{m-1}(\zeta_0, \dots, \zeta_{r-1})$.

*Induction:* we assume the induction hypotheses for some $j \in [m-1: r]$ and we want to prove the result for $j-1$. By equation (25), we have $\tilde{\mathfrak{U}}_j(z) = \mathcal{P}_j\tilde{\mathfrak{U}}_{j-1}(z)$ where

$$ \mathcal{P}_j = \begin{pmatrix} 1 & & & \\ & \ddots & & \\ & & \ddots & \\ & & & 1 \\ (0) & -\zeta_0 & \cdots & -\zeta_{r-1} \end{pmatrix} \in \mathcal{M}_{j+1,j}(\mathbb{C}). $$

We define $\mathcal{B}_{j-1}(\zeta_0, \dots, \zeta_{r-1}) = \mathcal{B}_j(\zeta_0, \dots, \zeta_{r-1})\mathcal{P}_j \in \mathcal{M}_{r,j}(\mathbb{C})$ then we have

$$ \mathcal{B}_{j-1}\tilde{\mathfrak{U}}_{j-1}(z) = \mathcal{B}_j\mathcal{P}_j\tilde{\mathfrak{U}}_{j-1}(z) = \mathcal{B}_j\tilde{\mathfrak{U}}_j(z) = \mathcal{B}\tilde{\mathfrak{U}}_{m-1}(z). $$

*Conclusion:* we define $\tilde{\mathcal{B}}$ by $\mathcal{B}_{r-1}$.

The function $\tilde{\mathcal{B}}$ is easily computable because $(\mathcal{P}_j)_j$ are just matrices of Gaussian elimination. $\square$

By (26), the intrinsic Kreiss-Lopatinskii determinant can be written

$$ \Delta(z) = \frac{\det(zK_{0,r-1}(z) - \mathcal{B}K_{0,m-1}(z))}{\det K_{0,r-1}(z)} = \frac{\det(zK_{0,r-1}(z) - \tilde{\mathcal{B}}K_{0,r-1}(z))}{\det K_{0,r-1}(z)} = \det(zI_r - \tilde{\mathcal{B}}). \quad (27) $$

The matrix $\tilde{\mathcal{B}}$ from Proposition 21 depends on coefficients $(\sigma_j(z))_j$. By (26), we have

$$ \tilde{\mathcal{B}}(\sigma_0(z), \dots, \sigma_{r-1}(z)) = \mathcal{B}K_{0,m-1}(z)K_{0,r-1}^{-1}(z). $$

Using any computer algebra system, we can compute the matrix $\tilde{\mathcal{B}}(\zeta_0, \dots, \zeta_{r-1})$ from the matrix $\mathcal{B}$ and then compute the coefficients $(\sigma_j(z))_j$ and replace $\zeta_j$ by $\sigma_j(z)$ for all $j \in [0: r-1]$. It provides that the computation of $\tilde{\mathcal{B}}(\zeta_0, \dots, \zeta_{r-1})$ can be done only once and then apply for different $z$ (and so different $(\sigma_j(z))_j$).

With (27), we again find the holomorphic property of $\Delta$ by Lemma 16 which states that $z \mapsto \sigma_j(z)$ is holomorphic on $\mathcal{U}$ for all $j \in [0: r-1]$.
but

$$ K_{0,3}(z_n) = \begin{pmatrix} 1 & 1 \\ \kappa_1(z_n) & \kappa_2(z_n) \\ \kappa_1^2(z_n) & \kappa_2^2(z_n) \\ \kappa_1^3(z_n) & \kappa_2^3(z_n) \end{pmatrix} \underset{n \to \infty}{\longrightarrow} K_{0,3}(z_\infty) = \begin{pmatrix} 1 & 0 \\ \kappa_1(z_\infty) & \kappa_1(z_\infty) \\ \kappa_1^2(z_\infty) & 2\kappa_1^2(z_\infty) \\ \kappa_1^3(z_\infty) & 3\kappa_1^3(z_\infty) \end{pmatrix}. $$

Consequently, the considered basis (2.4) of $\mathcal{E}^s(z)$ does not generally define a continuous mapping with respect to $z$.

In spite of the difficulty enlightened above, it turns out that $\mathcal{E}^s(z)$ is a continuous and even holomorphic vector bundle over $\mathcal{U}$ as it is discussed in [7, Thm. 4.3]. It is also proved that this vector bundle $\mathcal{E}^s(z)$ can even be continuously extended over $\bar{\mathcal{U}}$, thus considering $z \in \mathcal{S}$ as well (see also [37] for a similar property for the hyperbolic-parabolic PDE case). The main point therein is that for some $z_0 \in \mathcal{S}$, there may exists one (or several) root $\kappa_0(z_0)$ of (2.3) on $\mathcal{S}$. At such points $z_0$ the Lemma 2.1 (Hersh) does not hold anymore. This situation is depicted on the third line of Figure 2.1. For $z$ on $\mathcal{S}$, the space $\mathcal{E}^s(z)$ still is of dimension $r$ and we extend the notation $K_{i,j}(z)$. We can summarize the above discussion in the following theorem.

**THEOREM 2.6 ([7]).** — Under Assumption 1.2, the space $\mathcal{E}^s(z)$ is a holomorphic vector bundle over $\mathcal{U}$ and can be extended in a unique way to a continuous vector bundle over $\bar{\mathcal{U}}$.

**Remark 2.7.** — For the extension, the first difficulty is to select the roots of (2.3) coming from the inside, indeed, if there is a root on $\mathcal{S}$, it can be coming from the inside of $\mathbb{D}$, the outside or both (in case of multiplicity). In Section 4.2, we will explain the numerical strategy to select the good ones. The second difficulty is to prove the continuity of $\mathcal{E}^s(z)$ after the extension, it follows from the existence of a K-symmetrizer and is obtained e.g. in [7, Thm. 4.3]. As previously observed, $K_{i,j}(z)$ is generally not continuous with respect to $z$.

## 2.2. Intrinsic Kreiss–Lopatinskii determinant

In this section, we define properly formulas (1.10) and (1.11). Let us consider the $\mathcal{Z}$-transformed version of the boundary condition (1.7), that is

$$ z \begin{pmatrix} \tilde{U}_0(z) \\ \vdots \\ \tilde{U}_{r-1}(z) \end{pmatrix} - \mathcal{B} \begin{pmatrix} \tilde{U}_0(z) \\ \vdots \\ \tilde{U}_{m-1}(z) \end{pmatrix} = \begin{pmatrix} a_{-r} & \cdots & a_{-1} \\ \vdots & \ddots & \vdots \\ 0 & \cdots & a_{-r} \end{pmatrix} \begin{pmatrix} \tilde{g}_{-r}(z) \\ \vdots \\ \tilde{g}_{-1}(z) \end{pmatrix}. \quad (2.5) $$

Injecting the solution $(\tilde{U}_j(z))_{j \in \mathbb{N}} \in \mathcal{E}^s(z)$ to (2.2) into (2.5), we obtain a system of $r$ equations with $r$ scalar unknowns: they are the coefficients of $(\tilde{U}_j(z))_{j \in \mathbb{N}}$ written in the basis (2.4) of $\mathcal{E}^s(z)$.
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Illustration of Lemma 22</figcaption></figure>

4.2. **Computation of $\Delta(\mathbb{S})$.** Let us fix a $z_0 \in \mathbb{S}$. To compute $(\sigma_j(z_0))_j$, we need the $r$ roots $(\kappa_j(z_0))_j$ that come from the inside of the unit disk, see Remark 8. By the continuity of the roots of polynomial $P_{z_0}$ defined in (20) with respect to the parameter $z_0$, for each $\kappa_0(z_0)$ of multiplicity $\beta$ on the unit circle, for a sufficiently small $\varepsilon > 0$, there exists $\eta > 0$ such that for all $z \in B(z_0, \eta)$, the polynomial $P_z$ has exactly $\beta$ roots with multiplicity in $B(\kappa_0(z_0), \varepsilon)$. The explicit value of $\eta$ is given in the following statement.

**Lemma 22.** Let $z_0$ be on the unit circle. Let $\kappa_0(z_0) \in \mathbb{S}$ be a root of multiplicity $\beta$ of the polynomial $P_{z_0}$ defined in (20). Let $\varepsilon > 0$ be such that $\kappa_0(z_0)$ is the only root of $P_{z_0}$ in $\overline{B(\kappa_0(z_0), \varepsilon)}$ and set

$$ \eta = (1 + \varepsilon)^{-r} \min_{\kappa \in \partial B(\kappa_0(z_0), \varepsilon)} |P_{z_0}(\kappa)|. $$

Then for all $z \in B(z_0, \eta)$, the polynomial $P_z$ has exactly $\beta$ roots with multiplicity in $B(\kappa_0(z_0), \varepsilon)$.

The proof of Lemma 22 is a consequence of Rouché's theorem, comparing the number of zeros between $P_{z_0}$ and $P_z$ for $z$ close to $z_0$ in $B(\kappa_0(z_0), \varepsilon)$, the details of the proof are not given here, but we refer to [24] for a proof of Rouché's theorem. This lemma is illustrated in Figure 2 where $\kappa_0(z_0)$ is of multiplicity 3. The black points are related to $z_0$ and $(\kappa_j(z_0))_j$, and the gray points are related to $z$ and $(\kappa_j(z))_j$.

From the numerical point of view, for a multiple root $\kappa_0(z_0)$, one can take the smallest distance between two roots of $P_{z_0}$ as $\varepsilon$, take $z = (1 + \frac{\eta}{2})z_0 \in \mathcal{U}$. The value of $\eta$ is obtained discretizing the circle of radius $\varepsilon$ centered in $\kappa_0(z_0)$. By Lemma 3 (Hersh), there is no roots of $P_z$ on the unit circle, then one can count the roots in $B(\kappa_0(z_0), \varepsilon) \cap \mathbb{D}$ and $B(\kappa_0(z_0), \varepsilon) \cap \mathcal{U}$ to know the number of roots linked to $\kappa_0(z)$ that come from the inside and the outside of the unit disk. After selecting the roots $(\kappa_j(z))_j$ that come from the inside of the unit disk, one may compute their symmetric functions $(\sigma_j(z))_j$. By replacing the formal variables $(\zeta_j)_j$ of $\tilde{B}(\zeta_0, \dots, \zeta_{r-1})$ with $(\sigma_j(z))_j$, one may compute $\Delta(z)$ with expression (27). Instead of computing $(\kappa_j(z))_j$ for each $z$ on the unit circle independently, one may use the continuity of $(\kappa_j(z))_j$ with respect to $z$ in order to describe the movement of the roots $(\kappa_j(z))_j$ for $z \in \mathbb{S}$. After drawing the Kreiss-Lopatinskii curve, the winding number has to be computed in order to use Method 15. To do so, we use the geometric algorithm proposed by García Zapata and Díaz Martín in [13] and [14].

4.3. **Boundary condition: reconstruction procedure.** To define the boundary condition, we use the reconstruction procedure explained in [10]. The framework is the advection equation with a
**Remark 2.8.** — For $r = 2$ and a given value of $z$ (we skip for convenience the dependence in $z$ hereafter), if $\kappa_1 \neq \kappa_2$ so that the solution to (2.5) has the form $\alpha_1\kappa_1^j + \alpha_2\kappa_2^j$, then that solution is constrained by the system (2.5). The matricial form of that system reads

$$ \left( z \begin{pmatrix} 1 & 1 \\ \kappa_1 & \kappa_2 \end{pmatrix} - \mathcal{B} \begin{pmatrix} 1 & 1 \\ \kappa_1 & \kappa_2 \\ \vdots & \vdots \\ \kappa_1^{m-1} & \kappa_2^{m-1} \end{pmatrix} \right) \begin{pmatrix} \alpha_1 \\ \alpha_2 \end{pmatrix} = \begin{pmatrix} a_{-1} \tilde{g}_{-1} + a_{-2} \tilde{g}_{-2} \\ a_{-2} \tilde{g}_{-1} \end{pmatrix}. $$

The injectivity, whence invertibility, of the boundary condition is thus directly related to the property $\det(zK_{0,1} - BK_{0,m-1}(z)) \neq 0$, where $zK_{0,1} - BK_{0,m-1}(z) \in \mathcal{M}_{2,2}(\mathbb{C})$.

**DEFINITION 2.9** (Kreiss–Lopatinskii determinant). — The Kreiss–Lopatinskii determinant is the complex-valued function defined for $|z| \ge 1$ by:

$$ \Delta_{\mathrm{KL}}(z) = \det(zK_{0,r-1}(z) - BK_{0,m-1}(z)). $$

Despite the fact that the space $\mathcal{E}^s(z)$ is a holomorphic vector bundle over $\mathcal{U}$ and continuous over $\bar{\mathcal{U}}$ (Theorem 2.6), this determinant $\Delta_{\mathrm{KL}}$ is not holomorphic on $\mathcal{U}$. To retrieve those properties, we define the intrinsic Kreiss–Lopatinskii determinant $\Delta$ that we can motivate by the following informal discussion. The above Kreiss–Lopatinskii determinant is actually not well defined until we order in some way the roots $(\kappa_j(z))_{j=1,...,r}$ of (2.3). There are two points to emphasize. The first one is related to crossing roots and already discussed after Remark 2.5. The second one is that, outside crossing cases, being given any choice for the ordering of the roots (and thus of the vectors of the basis (2.4) for the vector bundle), there is in general no chance to obtain a holomorphicity property for the components of the matrix $K_{0,m-1}(z)$ over $\mathcal{U}$. For example, even the roots of $X^2 - z$ are not holomorphic w.r.t. $z \in \mathcal{U}$ because of the logarithm determination. On the other side, any symmetric functions of the roots $(\kappa_j(z))_{j=1,...,r}$ however are holomorphic because they can be obtained directly in terms of the coefficients of the polynomial (2.3). So except for crossing roots, the same holds for the quantity $\Delta_{\mathrm{KL}}(z)$ since the matrices $B$ and $\tilde{B}$ are constants and the determinant itself is a symmetric function.

A very natural way to reach the holomorphic property and go beyond the last difficulties consists in dividing $\Delta_{\mathrm{KL}}$ by the quantity $\det K_{0,r-1}(z)$. Hence, the same permutation or combination of the vectors of the basis (2.4) is involved in both computations. This intrinsic Kreiss–Lopatinskii determinant has already been introduced and studied in a particular case in [5].
misalignment between the space boundary and the discrete grid points

$$ egin{cases} \partial_t u + a \partial_x u = 0, & t \ge 0, x \in [x_\sigma, 1], \\ u(t, x_\sigma) = g(t), & t \ge 0, \\ u(0, x) = f(x), & x \in [x_\sigma, 1]. \end{cases} \tag{28} $$

Without loss of generality, we can assume that $x_\sigma = \sigma \Delta x$ with $\sigma \in [-\frac{1}{2}, \frac{1}{2}]$ (as it is explained in [4]). Let us introduce $x_j$ for $j \Delta x$ and $t^n$ for $n \Delta t$ when $j \ge -r$ and $n \ge 0$. Let $n \in \mathbb{N}$ be a fixed time. The solution $u$ of (28) (assumed here to be smooth enough) satisfies

$$ \frac{1}{\Delta x} \int_{x_j - \frac{\Delta x}{2}}^{x_j + \frac{\Delta x}{2}} u(t^n, y) dy = \frac{1}{\Delta x} \int_{x_j - \frac{\Delta x}{2}}^{x_j + \frac{\Delta x}{2}} \sum_{k=0}^{d-1} \partial_x^k u(t^n, x_\sigma) \frac{(y - x_\sigma)^k}{k!} dy + O(\Delta x^d) \tag{29} $$

using a Taylor expansion of order $d$. Now, let us take a solution $(U_j^n)_{j \ge 0}$ of a scheme of the form (2) approximating $u$. Using (29), we want to define the $r$ ghost points $(U_j^{n+1})_{-r \le j \le -1}$. The approximation of equation (29) reads, for $j \ge -r$,

$$ U_j^n \approx \sum_{k=0}^{d-1} \partial_x^k u(t^n, x_\sigma) \left( \frac{(j + \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} - \frac{(j - \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} \right). \tag{30} $$

On the one hand, we use the PDE to convert space derivatives into time derivatives until an index $k_d < d$. The index $k_d$ allows us to know only the first derivatives of the boundary datum $g$ and use extrapolation for the rest. For the advection equation, we have, for all $k \le k_d$,

$$ \partial_x^k u(t^n, x_\sigma) = (-a)^{-k} \partial_t^k u(t^n, x_\sigma) = (-a)^{-k} g^{(k)}(t^n). $$

Equation (30) becomes, for $j \ge -r$,

$$ U_j^n \approx \sum_{k=0}^{k_d} (-a)^{-k} g^{(k)}(t^n) \left( \frac{(j + \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} - \frac{(j - \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} \right) + \sum_{k=k_d+1}^{d-1} \partial_x^k u(t^n, x_\sigma) \left( \frac{(j + \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} - \frac{(j - \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} \right). \tag{31} $$

On the other hand, we need to define $(\partial_x^k u(t^n, x_\sigma))_{k=k_d+1}^{d-1}$, but using (31) for $j \in [0 : d - k_d - 2]$, we can deduce the unknowns $(U_j^n)_{-r \le j \le -1}$. Writing $\mathcal{U}_-$ = $(U_{-r}^n, \dots, U_{-1}^n)^T$, $\mathcal{U}_+$ = $(U_0^n, \dots, U_{d-k_d-2}^n)^T$ and $\Theta^n = (\partial_x^{k_d+1} u(t^n, x_\sigma), \dots, \partial_x^{d-1} u(t^n, x_\sigma))^T$, we have a condensed formulation of (31):

$$ \left\{ \begin{array}{l} \mathcal{U}_{-} = S_{-}^{n} + \mathcal{Y}_{-}^{\Theta^{n}}, \\ \mathcal{U}_{+} = S_{+}^{n} + \mathcal{Y}_{+}^{\Theta^{n}}, \end{array} \right. \tag{32} $$

where $S_-^n \in \mathbb{R}^r$, $S_+^n \in \mathbb{R}^{d-k_d-1}$, $\mathcal{Y}_- \in \mathcal{M}_{r,d-k_d-1}(\mathbb{R})$ and $\mathcal{Y}_+ \in \mathcal{M}_{d-k_d-1}(\mathbb{R})$ with

$$ \left\{ \begin{array}{ll} (S_-^n)_i = \sum_{k=0}^{k_d} (-a)^{-k} g^{(k)}(t^n) \left( \frac{(-i + \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} - \frac{(-i - \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} \right) & \text{for } i \in [1:r] \\ (S_+^n)_i = \sum_{k=0}^{k_d} (-a)^{-k} g^{(k)}(t^n) \left( \frac{(i - 1 + \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} - \frac{(i - 1 - \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} \right) & \text{for } i \in [1:d - k_d - 1] \\ (\mathcal{Y}_-)_i, j = \left( \frac{(- (r - i + 1) + \frac{1}{2} - \sigma)^{j+k_d+1}}{(j+k_d+1)!} - \frac{(- (r - i + 1) - \frac{1}{2} - \sigma)^{j+k_d+1}}{(j+k_d+1)!} \right) & \text{for } i \in [1:r], \\ (\mathcal{Y}_+)_i, j = \left( \frac{(i - 1 + \frac{1}{2} - \sigma)^{j+k_d+1}}{(j+k_d+1)!} - \frac{(i - 1 - \frac{1}{2} - \sigma)^{j+k_d+1}}{(j+k_d+1)!} \right) & \text{for } i, j \in [1:d - k_d - 1] \end{array} \right. $$
**DEFINITION 2.10** (Intrinsic Kreiss–Lopatinskii determinant). — The intrinsic Kreiss–Lopatinskii determinant is the complex-valued function defined for $|z| \ge 1$ by:

$$ \Delta(z) = \frac{\Delta_{\text{KL}}(z)}{\det K_{0,r-1}(z)}. \qquad (2.6) $$

Let us note that the intrinsic Kreiss–Lopatinskii determinant can be rewritten

$$ \begin{aligned} \Delta(z) &= \frac{\det(zK_{0,r-1}(z) - BK_{0,m-1}(z))}{\det K_{0,r-1}(z)} \\ &= z^r \det \left( I_r - \frac{BK_{0,m-1}(z)K_{0,r-1}(z)^{-1}}{z} \right). \end{aligned} \qquad (2.7) $$

To conclude with these definitions, let us state a little more about the *Uniform Kreiss–Lopatinskii Condition*. With the above notations and additionally to the invertibility of $zK_{0,r-1}(z) - BK_{0,m-1}(z)$, it corresponds to the existence of a constant $C > 0$ such that for any $z \in \bar{U}$, any $\tilde{U} \in \mathcal{E}^s(z)$ solution to (2.5) satisfies the uniform estimate $\|\tilde{U}\| \le C\|\tilde{g}\|$. From the Parseval identity for the $\mathcal{Z}$-transform, this inequality directly gives the first necessary half-part of the strong stability estimate (1.9). We refer the reader to [22] for a more detailed presentation.

## 2.3. Main results

Theorem 2.11 is our main theoretical result. It states that the intrinsic Kreiss–Lopatinskii determinant has the same regularity properties as $\mathcal{E}^s(z)$, see Theorem 2.6.

**THEOREM 2.11** (Smoothness of the intrinsic Kreiss–Lopatinskii determinant). — Assume 1.2. The intrinsic Kreiss–Lopatinskii determinant $\Delta$ is holomorphic on $U$ and continuous on $\bar{U}$.

By (2.6), the function $\Delta$ shares the same zeros with the Kreiss–Lopatinskii determinant $\Delta_{\text{KL}}$, so that it can be used as an alternative in the Uniform Kreiss–Lopatinskii Condition, see Theorem 1.4 (Kreiss). Another property, important for the forthcoming applications, lies in the next Corollary 2.13 and involves the following important geometrical object:

**DEFINITION 2.12.** — The Kreiss–Lopatinskii curve $\Delta(\mathcal{S})$ is the closed complex parameterized curve

$$ \Delta(\mathcal{S}) = \{\theta \in [0, 2\pi] \longmapsto \Delta(e^{i\theta})\}. $$
Eliminating the space derivatives $\Theta^n$ in (32) gives us the following boundary condition:

$$ \mathcal{U}_- = \mathcal{Y}_- \mathcal{Y}_+^{-1} \mathcal{U}_+ + \mathcal{S}_-^n - \mathcal{Y}_- \mathcal{Y}_+^{-1} \mathcal{S}_+^n. \quad (33) $$

Equation (33) is exactly the boundary equations (5) of the scheme which define the r ghost points of the scheme. To write the boundary condition as equation (7) with expression (8), we identify

$$ B \stackrel{\text{def}}{=} \mathcal{Y}_- \mathcal{Y}_+^{-1} \text{ and } \begin{pmatrix} g_{-r}^n \\ \vdots \\ g_{-1}^n \end{pmatrix} \stackrel{\text{def}}{=} \mathcal{S}_-^n - \mathcal{Y}_- \mathcal{Y}_+^{-1} \mathcal{S}_+^n. $$

As in [10], $\mathcal{R}^{d,k_d}$ denotes the reconstruction procedure where $d$ is the order of consistency of the method and $k_d$ the index when we change from time derivatives to extrapolation. For example, the reconstruction procedure $\mathcal{R}^{3,0}$ for $r = 2$ and $\sigma = 0.4$ leads to

$$ \mathcal{Y}_- = \begin{pmatrix} -\frac{12}{5} & \frac{1753}{600} \\ -\frac{7}{5} & \frac{613}{600} \end{pmatrix}, \mathcal{Y}_+ = \begin{pmatrix} -\frac{2}{5} & \frac{73}{600} \\ \frac{3}{5} & \frac{133}{600} \end{pmatrix} \text{ and } B = \begin{pmatrix} \frac{1371}{97} & \frac{526}{97} \\ \frac{554}{97} & \frac{143}{97} \end{pmatrix} $$

4.4. **Example of O3 scheme.** As it is done in [10], we want to find the stability area for the O3 scheme defined, for $j \in \mathbb{N}$ and $n \in \mathbb{N}$, by

$$ U_{j+1}^{n+1} = \left(\frac{\lambda^3}{6} - \frac{\lambda}{6}\right) U_{j-2}^{n} + \left(\lambda + \frac{\lambda^2}{2} - \frac{\lambda^3}{2}\right) U_{j-1}^{n} + \left(1 - \frac{\lambda}{2} - \lambda^2 + \frac{\lambda^3}{2}\right) U_{j}^{n} + \left(\frac{\lambda^2}{2} - \frac{\lambda^3}{6} - \frac{\lambda}{3}\right) U_{j+1}^{n} \quad (34) $$

The O3 scheme is a scheme with $r = 2$ and $p = 1$ and is Cauchy-stable for $\lambda \in ]0, 1]$. The reconstruction $\mathcal{R}^{3,0}$ for the O3 scheme and $\sigma = 0.4$ leads to

$$ B = \begin{pmatrix} \frac{1371}{97} & \frac{526}{97} & 0 \\ \frac{554}{97} & \frac{143}{97} & 0 \end{pmatrix} \text{ and } \mathcal{B} = \begin{pmatrix} \frac{180\lambda^2}{97} + \frac{277\lambda}{97} + 1 & \frac{120\lambda^2}{97} + \frac{23\lambda}{97} & 0 \\ \frac{263\lambda^3}{582} + \frac{\lambda^2}{2} + \frac{14\lambda}{291} & \frac{217\lambda^3}{291} - \lambda^2 - \frac{217\lambda}{291} + 1 & -\frac{\lambda^3}{6} + \frac{\lambda^2}{2} - \frac{\lambda}{3} \end{pmatrix}. $$

Using the reformulation (27) of the Kreiss-Lopatinskii determinant, we have

$$ \tilde{\mathcal{B}}(\zeta_0, \zeta_1) = \mathcal{B} \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ -\zeta_0 & -\zeta_1 \end{pmatrix} = \begin{pmatrix} \frac{180\lambda^2}{97} + \frac{277\lambda}{97} + 1 & \frac{120\lambda^2}{97} + \frac{23\lambda}{97} & 0 \\ \frac{(263+97\zeta_0)\lambda^3}{582} + \frac{(1-\zeta_0)\lambda^2}{2} + \frac{(14+97\zeta_0)\lambda}{291} & \frac{(434+97\zeta_1)\lambda^3}{582} - \frac{\lambda^2}{2} - \frac{(217-97\zeta_1)\lambda}{291} + 1 & 0 \end{pmatrix}. $$

For example, for $\sigma = 0.4$, Figure 3 shows that the O3 scheme with $\mathcal{R}^{3,0}$ boundary is stable for $\lambda = 0.4$ (because $r - \text{Ind}_{\Delta(\mathcal{S})}(0) = 0$) and is unstable for $\lambda = 0.9$ (because $r - \text{Ind}_{\Delta(\mathcal{S})}(0) = 1$).

We can draw the same figure as the Figure 4 of [10] but instead of using a computation of the spectral radius of the truncated quasi-Toeplitz matrix, we use our strategy of counting the number of instability modes, see Figure 4 which is much more reliable (since it is parameter free) and efficient. In Figure 4, every area stamped with 0 is a domain where the O3 scheme is stable. The odd pattern for very small $\lambda$ (approximatively between 0 and 0.01) of Figure 4 may be due to difficulties for computing the winding number. Indeed, for very small values of $\lambda$, the Kreiss-Lopatinskii determinant is really close to the origin and even with a refinement (for more details on this procedure see the next example), the computation of the winding number may become inaccurate, which is not a problem in practice since it would correspond to very small, then unusable, time steps.
Using the residue theorem<sup>(1)</sup> thanks to Theorem 2.11, we obtain the following result.

**COROLLARY 2.13** (Number of zeros of the intrinsic Kreiss–Lopatinskii determinant). — Assume 1.2. If $0 \notin \Delta(\mathbb{S})$ then the equation $\Delta(z) = 0$ has exactly $r - \text{Ind}_{\Delta(\mathbb{S})}(0)$ zeros in $\mathcal{U}$.

Here above and in all the paper, $\text{Ind}_{\Delta(\mathbb{S})}(0)$ denotes the winding number of the origin with respect to the closed oriented curve $\Delta(\mathbb{S})$ (see [29] for a definition of the winding number). This corollary helps us to establish an efficient and practical method to study the stability of a given IBVP through Theorem 1.4 (Kreiss). In particular, the low computational cost of the following procedure is very appealing for the study of parameterised IBVP's, see Section 4.

**METHOD 2.14** (Uniform Kreiss–Lopatinskii Condition check). — There are two different cases:

* if $0 \in \Delta(\mathbb{S})$, then there exists $z_0 \in \mathbb{S}$ such that $\Delta(z_0) = 0$.

* if $0 \notin \Delta(\mathbb{S})$, $\Delta$ does not vanish on $\mathbb{S}$ and it has $r - \text{Ind}_{\Delta(\mathbb{S})}(0)$ zeros in $\mathcal{U}$ by Corollary 2.13. It follows that if $\text{Ind}_{\Delta(\mathbb{S})}(0) = r$ then the scheme is stable. Otherwise the scheme is unstable.

In summary, by Theorem 1.4 (Kreiss) and since Uniform Kreiss–Lopatinskii Condition is fulfilled if and only if the Kreiss–Lopatinskii determinant has no zero in $\bar{\mathcal{U}}$, Method 2.14 can be used to conclude that the scheme is stable or not. Some illustrations for the O3 scheme and the fifth-order Lax–Wendroff scheme follow in Section 4.

### 3. Proof of Theorem 2.11 and Corollary 2.13

#### 3.1. Constant-recursive sequence of order $r$

For each $z \in \mathcal{U}$, we denote $P_z$ the polynomial linked to the characteristic equation (2.3), i.e.

$$ P_z(\kappa) = a_p\kappa^{r+p} + \dots + a_1\kappa^{r+1} + (a_0 - z)\kappa^r + a_{-1}\kappa^{r-1} + \dots + a_{-r+1}\kappa + a_{-r}. \quad (3.1) $$

By Lemma 2.1 (Hersh), the polynomial $P_z(\kappa)$ can be factorized into two polynomials: one with the $r$ roots in $\mathbb{D}$, denoted $R_z(\kappa)$ and one with the $p$

<footnote><sup>(1)</sup> All the complex analysis results can be found in [29].</footnote>
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Curve Δ(S) for O3 scheme for σ = 0.4, for λ ∈ {0.4, 0.9} with reconstruction boundary R$^{3,0}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Number of zeros of the Kreiss-Lopatinskii determinant of O3 scheme with different reconstruction boundaries for λ ∈ ]0, 1] and σ ∈ ] − 0 . 5 , 0 . 5[.</figcaption></figure>
roots in $\mathcal{U}$, denoted $Q_z(\kappa)$. We know that the coefficients of $P_z$ are holomorphic in $z$. We already said that the basis (2.4) is not holomorphic because the roots $\kappa$ are not. In the next result, we prove that the symmetric functions of the $r$ roots $\kappa$ living in $\mathbb{D}$ are indeed holomorphic in $\mathcal{U}$, in other words, the coefficients of $R_z$ are holomorphic in $\mathcal{U}$.

**LEMMA 3.1.** — For all $z \in \mathcal{U}$, the polynomial $R_z(X) = \prod_{j=1}^r (X - \kappa_j(z))$ has holomorphic coefficients in $\mathcal{U}$, where $(\kappa_j(z))_{j=1}^r$ are the $r$ roots (with multiplicity) in $\mathbb{D}$ of (2.3).

*Proof.* — We use the Dunford–Taylor formula with $C_z$ the companion matrix of the polynomial (2.3):

$$ \Pi(z) = \frac{1}{2\pi} \int_{\mathcal{S}} (\zeta I_{r+p} - C_z)^{-1} d\zeta $$

It is the projection along $\mathbb{E}_s(z) = \ker \prod_{j=1}^r (C_z - \kappa_j(z))$ onto $\mathbb{E}_u(z) = \ker \prod_{j=r+1}^{r+p} (C_z - \kappa_j(z))$ where $(\kappa_j(z))_{j=r+1}^{r+p}$ are the roots of (2.3) in $\mathcal{U}$, because $(\kappa_j(z))_{j=1}^r$ are surrounded by $\mathcal{S}$ and $(\kappa_j(z))_{j=r+1}^{r+p}$ are not. The projector $\Pi(z)$ is holomorphic on $\mathcal{U}$ since it is a holomorphic parameter integral. We have $C_z \circ \Pi(z)|_{\mathbb{E}_u(z)} = 0$ and $C_z \circ \Pi(z)|_{\mathbb{E}_s(z)} = C_z$, then the characteristic polynomial of $C_z \circ \Pi(z)$ is $X^p R_z(X)$ because $\mathbb{C}^{r+p} = \mathbb{E}_s(z) \oplus \mathbb{E}_u(z)$. The function $z \mapsto C_z \circ \Pi(z)$ is holomorphic on $\mathcal{U}$, then the coefficients of its characteristic polynomial are too. It concludes the proof. $\square$

## 3.2. Hermite interpolation

To prove the holomorphic properties of $\Delta$, by (2.7), it is sufficient to study the function $z \mapsto K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$. To simplify this study, for $z \in \bar{\mathcal{U}}$, we introduce on the linear map

$$ \varphi_z : Q \in \mathbb{C}_{m-1}[X] \mapsto \varphi_z(Q) \in \mathbb{C}_{r-1}[X] \quad (3.2) $$

where $\varphi_z(Q)$ is the Hermite interpolation polynomial of degree less than $r-1$ defined by the value

$$ (Q(\kappa_1(z)), \dots, Q^{(\beta_1-1)}(\kappa_1(z)), Q(\kappa_2(z)), \\ \dots, Q^{(\beta_2-1)}(\kappa_2(z)), \dots, Q^{(\beta_M-1)}(\kappa_M(z))), $$

where the $\kappa$'s are the same as in (2.4). The link between the product of matrices $K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$ and $\varphi_z$ is given by:

**LEMMA 3.2.** — For all $z \in \bar{\mathcal{U}}$, the matrix $(K_{0,m-1}(z)K_{0,r-1}^{-1}(z))^T$ is the representation in the canonical basis of $\varphi_z$ defined in (3.2).
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Number of zeros of the Kreiss-Lopatinskii determinant of LW5 scheme with different reconstruction boundaries for $\lambda \in ]0, 1]$ and $\sigma = 0.4$.</figcaption></figure>

4.5. **Example of Lax-Wendroff 5.** The fifth-order Lax-Wendroff scheme, called LW5, has been proposed in [29] and can be written, for all $n \in \mathbb{N}$, for all $j \in \mathbb{N}$, as

$$ U_j^{n+1} = \frac{\lambda(\lambda-2)(\lambda-1)(\lambda+1)(\lambda+2)}{120} U_{j-3}^n - \frac{\lambda(\lambda-1)(\lambda-3)(\lambda+1)(\lambda+2)}{24} U_{j-2}^n + \frac{\lambda(\lambda-2)(\lambda-3)(\lambda+1)(\lambda+2)}{12} U_{j-1}^n \\ + \left(1 - \frac{\lambda(\lambda^4-3\lambda^3-5\lambda^2+15\lambda+4)}{12}\right) U_j^n + \frac{\lambda(\lambda-1)(\lambda-2)(\lambda-3)(\lambda+2)}{24} U_{j+1}^n - \frac{\lambda(\lambda-1)(\lambda-2)(\lambda-3)(\lambda+1)}{120} U_{j+2}^n. \quad (35) $$

This scheme LW5 is Cauchy-stable for $\lambda \in ]0, 1]$.

Figure 5 illustrates the computation of the number of instabilities for LW5 for different reconstruction boundaries where $\sigma = 0.4$ with respect to $\lambda \in ]0, 1]$. As in the previous example, it may happen that the Kreiss-Lopatinskii curve is too close to the origin, the winding number of the origin cannot be computed correctly. Following the geometric algorithm proposed by García Zapata in [14], a refinement of the discretization then improves the effective computation of the winding number. Figure 6 represents such refinement with close-up close to the origin. However, even with this strategy, for very small values of $\lambda$, we cannot refine more than the machine's precision, that is why there is still some odd pattern for very small $\lambda$ in Figure 4, Figure 5 and Figure 7. As we already discussed, such very small time step are however not used in practice. The stability area with respect to both parameters $\lambda$ and $\sigma$ are drawn in Figure 7, considering again successively various reconstruction boundary conditions.

All the figures can be easily computed in Python with the common NumPy [20] library and the SymPy [31] library for the computer algebra system. The algorithm is really efficient. For each subfigure of Figure 4, the 1600 runs takes less than a couple of minutes of computation achieved on a standard laptop. Moreover, our procedure provides sharp results. In particular, contrary to numerical investigations of stability, which are based on the computation of the spectral radius, no arbitrary truncation of (quasi-)Toeplitz matrices is needed.
Proof. — The Hermite interpolation make appear the following matrix

$$
H_{0,j}(z) = \underbrace{\begin{pmatrix}
1 & 0 & 0 & \cdots & 1 & 0 & \cdots & 1 & 0 & \cdots \\
\kappa_1 & 1 & 0 & \cdots & \kappa_2 & 1 & \cdots & \kappa_M & 1 & \cdots \\
\kappa_1^2 & 2\kappa_1 & 2 & \cdots & \kappa_2^2 & 2\kappa_2 & \cdots & \kappa_M^2 & 2\kappa_M & \cdots \\
\kappa_1^3 & 3\kappa_1^2 & 6\kappa_1 & \cdots & \kappa_2^3 & 3\kappa_2^2 & \cdots & \kappa_M^3 & 3\kappa_M^2 & \cdots \\
\kappa_1^4 & 4\kappa_1^3 & 12\kappa_1^2 & \cdots & \kappa_2^4 & 4\kappa_2^3 & \cdots & \kappa_M^4 & 4\kappa_M^3 & \cdots \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \\
\kappa_1^j & j\kappa_1^{j-1} & j(j-1)\kappa_1^{j-2} & \cdots & \kappa_2^j & j\kappa_2^{j-1} & \cdots & \kappa_M^j & j\kappa_M^{j-1} & \cdots
\end{pmatrix}}_{\text{\underline{\beta}_1 columns linked to } \kappa_1}
\underbrace{\begin{pmatrix}
\kappa_2 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\
\kappa_2^2 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\
\kappa_2^3 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\
\kappa_2^4 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \\
\kappa_2^j & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots
\end{pmatrix}}_{\text{\underline{\beta}_2 columns linked to } \kappa_2}
\underbrace{\begin{pmatrix}
\kappa_M & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\
\kappa_M^2 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\
\kappa_M^3 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\
\kappa_M^4 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \\
\kappa_M^j & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots
\end{pmatrix}}_{\text{\underline{\beta}_M columns linked to } \kappa_M}
$$

The representation of $\varphi_z$ in the canonical basis is $(H_{0,m-1}(z)H_{0,r-1}^{-1}(z))^T$.
Besides, there exists an invertible matrix $M(z) \in \mathcal{M}_r(\mathbb{C})$ such that $K_{0,j}(z) = H_{0,j}(z)M(z)$. Therefore, we have

$$
\begin{align*}
K_{0,m-1}(z) K_{0,r-1}^{-1}(z) &= H_{0,m-1}(z) M(z) M(z)^{-1} H_{0,r-1}^{-1}(z) \\
&= H_{0,m-1}(z) H_{0,r-1}^{-1}(z).
\end{align*}
$$

The result follows.

<figure><img src="image_6.png" /></figure>

PROPOSITION 3.3. — The function $z \mapsto \varphi_z$ is holomorphic on $\mathcal{U}$.

*Proof.* — For all $k \in [0:m-1]$, we want every coefficient of the polynomial $\varphi_z(X^k)$ to be holomorphic on $\mathcal{U}$. Writing $\varphi_z(X^k)(x) = \sum_{j=0}^{r-1} \alpha_{j,k}(z)x^j$, we know that

$$
\forall j \in [0 : r-1], \quad j! \alpha_{j,k}(z) = \partial_x^j \varphi_z(X^k)(x) \Big|_{x=0}. \quad (3.3)
$$

By the error of Hermite interpolation (see [25]), we have

$$
\varphi_z(X^k)(x) - x^k = \frac{1}{2i\pi} \int_{\mathbb{S}} \frac{\zeta^k R_z(x)}{(x - \zeta) R_z(\zeta)} d\zeta \quad (3.4)
$$

where $R_z(X)$ is defined in Lemma 3.1. Differentiating equation (3.4) (with
the Leibniz product rule), one obtains

$$
j! \alpha_{j,k}(z) = k! \delta_k^j + \sum_{s=0}^{j} \binom{j}{s} R_z^{(j-s)}(0) \frac{1}{2i\pi} \int_{\mathbb{S}} \frac{-s! \zeta^{k-s-1}}{R_z(\zeta)} d\zeta. \quad (3.5)
$$

By Lemma 3.1 and the holomorphicity of parameter-dependent integrals,
the function $z \mapsto \alpha_{j,k}(z)$ is holomorphic on $\mathcal{U}$ for all $j \in [0 : r-1]$ and
$k \in [0 : m-1]$. The proof is now complete.
$\square$

<figure><img src="image_7.png" /></figure>

PROPOSITION 3.4. — The function $z \mapsto \varphi_z$ is continuous on $\bar{\mathcal{U}}$.
<figure><img src="image_8.png" /><figcaption>FIGURE 6. (a) representation of $\Delta(\mathcal{S})$ for LW5 with the boundary condition $\mathcal{R}^{6,1}$, for $\lambda = 0.01$ and $\sigma = 0$, zoom (b) without refinement (c) and with refinement (d).</figcaption></figure>

REFERENCES

[1] A. Arnold, M. Ehrhardt, and I. Sofronov. Discrete transparent boundary conditions for the Schrödinger equation: fast calculation, approximation, and stability. *Communications in Mathematical Sciences*, 1(3):501–556, 2003.

[2] R. M. Beam and R. F. Warming. The asymptotic spectra of banded Toeplitz and quasi-Toeplitz matrices. *SIAM Journal on Scientific Computing*, 14(4):971–1006, 1993.

[3] S. Benzoni-Gavage and D. Serre. *Multi-dimensional hyperbolic partial differential equations: First-order Systems and Applications*. OUP Oxford, 2006.

[4] B. Boutin, P. Le Barbenchon, and N. Seguin. On the stability of totally upwind schemes for the hyperbolic initial boundary value problem. submitted, 2022.

[5] B. Boutin, T. Nguyen, A. Sylla, S. Tran-Tien, and J.-F. Coulombel. High order numerical schemes for transport equations on bounded domains. *ESAIM: Proceedings and Surveys*, 70:84–106, 2021.

[6] J.-F. Coulombel. Stability of finite difference schemes for hyperbolic initial boundary value problems. In *HCDTE lecture notes. Part I. Nonlinear hyperbolic PDEs, dispersive and transport equations*, volume 6 of *AIMS Ser. Appl. Math.*, page 146. Am. Inst. Math. Sci. (AIMS), Springfield, MO, 2013.

[7] J.-F. Coulombel. Transparent numerical boundary conditions for evolution equations: derivation and stability analysis. *Ann. Fac. Sci. Toulouse, Math. (6)*, 28(2):259–327, 2019.

[8] R. Courant, K. Friedrichs, and H. Lewy. Über die partiellen Differenzengleichungen der mathematischen Physik. *Mathematische Annalen*, 100(1):32–74, 1928.

[9] J. Crank and P. Nicolson. A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. *Mathematical Proceedings of the Cambridge Philosophical Society*, 43(1):50–67, 1947.

[10] G. Dakin, B. Després, and S. Jaouen. Inverse Lax-Wendroff boundary treatment for compressible Lagrange-remap hydrodynamics on cartesian grids. *Journal of Computational Physics*, 353:228–257, 2018.
*Proof.* — Because Lemma 2.1 (Hersh) does not hold anymore for $z \in \mathbb{S}$, the roots $\kappa(z)$ of characteristic equation (2.3) can be on the unit circle $\mathbb{S}$. To prove the continuity of $z \mapsto \alpha_{j,k}(z)$, we use equation (3.5) but replacing $\mathbb{S}$ by $\mathbb{S}_\varepsilon \stackrel{\text{def}}{=} \{z \in \mathbb{C}, |z| = 1 + \varepsilon\}$ for $\varepsilon > 0$. Using the continuity of parameter-dependent integrals and the continuity of the roots $\kappa(z)$ of characteristic equation (2.3), we obtain the continuity of the coefficients of $R_z$ and thus the function $z \mapsto \alpha_{j,k}(z)$ is continuous on $\bar{\mathcal{U}}$ for all $j \in [0 : r - 1]$ and $k \in [0 : m - 1]$. The proof is now complete. $\square$

**PROPOSITION 3.5.** — The function $z \mapsto \varphi_z$ is bounded on $\bar{\mathcal{U}}$.

*Proof.* — Equation (3.5) can give a bound of every components of $K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$. Indeed, using Rouché's theorem, as in the proof of Lemma 2.1 (Hersh), we can see that for $|z| > R$ for a certain $R$, all the roots $\kappa(z)$ of the characteristic equation (2.3) satisfy $|\kappa(z)| < \frac{1}{2}$. Then, for $|z| > R$, one can have

$$ \left| \frac{1}{2i\pi} \int_{\mathbb{S}} \frac{-s! \zeta^{k-s-1}}{R_z(\zeta)} d\zeta \right| = \left| \frac{s!}{2i\pi} \int_0^{2\pi} \frac{e^{i\theta(k-s-1)}}{\prod_{j=1}^r (e^{i\theta} - \kappa_j(z))} d\theta \right| \\ \le \frac{s!}{2\pi} \int_0^{2\pi} \frac{1}{|1 - \frac{1}{2}|^r} d\theta \le s! 2^r. $$

By Gauss-Lucas theorem, the roots of all the derivatives of $R_z$ are in $\mathbb{D}$ for all $|z| > 1$, it follows that $R_z^{(j-s)}(0)$ is bounded independently of $z$. Then for $|z| > R$, the quantity $K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$ is bounded. Moreover, by Proposition 3.4, the quantity $K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$ is bounded on the compact set $\{1 \le |z| \le R\}$. The proof is now complete. $\square$

### 3.3. Conclusion

*Proof of Theorem 2.11.* — By Lemma 3.2, the continuity and holomorphicity properties of $z \mapsto \varphi_z$ provided in Propositions 3.3 and 3.4 are shared by the function $z \mapsto K_{0,m-1}(z)K_{0,r-1}^{-1}(z)$. The expression of the intrinsic determinant (2.7) concludes the proof. $\square$

The next proof is close to the proof of the Corollary 15 of [5]. We reproduce it here for completeness.

*Proof of Corollary 2.13.* — Let us define the following function

$$ \tilde{\Delta} : z \in \mathbb{D}^* \mapsto \Delta(1/z) \in \mathbb{C}. $$

By Theorem 2.11, the function $\tilde{\Delta}$ is meromorphic on $\mathbb{D}$ with one only pole in 0 and is continuous on $\mathbb{D} \setminus \{0\}$. By Proposition 3.5, the function $z \mapsto$
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Number of zeros of the Kreiss-Lopatinskii determinant of LW5 scheme with different reconstruction boundaries for $\lambda \in ]0, 1]$ and $\sigma \in ] -0.5, 0.5[.$</figcaption></figure>
$K_{0,m-1}(1/z)K_{0,r-1}(1/z)^{-1}$ is bounded on $\bar{\mathbb{D}}$, it follows that 0 is a pole of order $r$ of the function $\tilde{\Delta}$. The residue theorem applied on $\tilde{\Delta}$ with the path $S$ gives the following equality:

$$ \mathrm{Ind}_{\tilde{\Delta}(S)}(0) = \#\mathrm{zeros}_{\tilde{\Delta}}(\bar{\mathbb{D}}) - \#\mathrm{poles}_{\tilde{\Delta}}(\bar{\mathbb{D}}). $$

It follows that

$$ \#\mathrm{zeros}_{\tilde{\Delta}}(\mathcal{U}) = r - \mathrm{Ind}_{\tilde{\Delta}(S)}(0). \quad \square $$

# 4. Numerical results

The results presented below can be reproduced from the code available in the GitHub repository [31], it can be used as a Python library, following the link: https://doi.org/10.5281/zenodo.7773742.

## 4.1. New formulation of $\tilde{\Delta}$

In [5], an explicit formula of the Kreiss-Lopatinskii determinant is given. Unfortunately to reduce the boundary matrix $\mathcal{B}$, the characteristic equation of degree $r+p$ is used to find a final matrix of size $r \times (r+p)$ and, since in the present analysis $p \neq 0$, the matrix is not square. To skirt that problem we will use the polynomial $R_z$ defined in Lemma 3.1 instead of using the complete characteristic equation (3.1). It reads also

$$ R_z(X) = \prod_{j=1}^{r} (X - \kappa_j(z)) = X^r + \sigma_{r-1}(z)X^{r-1} + \dots + \sigma_1(z)X + \sigma_0(z) $$

where $(\sigma_j(z))_j$ are the symmetric functions of $(\kappa_j(z))_j$.

Because $\tilde{U}_j(z)$ is in $\mathcal{E}^s(z)$ and can be expressed in the basis (2.4), we have, for all $j \in \mathbb{N}$,

$$ \tilde{U}_{j+r}(z) + \sigma_{r-1}(z)\tilde{U}_{j+r-1}(z) + \dots + \sigma_1(z)\tilde{U}_{j+1}(z) + \sigma_0(z)\tilde{U}_j(z) = 0. \quad (4.1) $$

**Notation 4.1.** — For all $j \in \mathbb{N}$, $\tilde{\mathfrak{U}}_j(z)$ is the vector $(\tilde{U}_0(z) \cdots \tilde{U}_j(z))^T$ of size $j+1$.

**PROPOSITION 4.2.** — Let $\mathcal{B} \in \mathcal{M}_{r,m}(\mathbb{C})$. There exists a function $\tilde{\mathcal{B}}: \mathbb{C}^r \to \mathcal{M}_{r,r}(\mathbb{C})$ constructible such that, for all $z \in \tilde{U}$, we have

$$ \mathcal{B}\tilde{\mathfrak{U}}_{m-1}(z) = \tilde{\mathcal{B}}(\sigma_0(z), \dots, \sigma_{r-1}(z))\tilde{\mathfrak{U}}_{r-1}(z) \quad (4.2) $$

where $(\tilde{U}_j(z))_j$ satisfies (4.1) for all $j \in \mathbb{N}$.
[11] M. Ehrhardt. Absorbing boundary conditions for hyperbolic systems. *Numer. Math.*, Theory Methods Appl., 3(3):295-337, 2010.

[12] B. Engquist and A. Majda. Absorbing boundary conditions for the numerical simulation of waves. *Mathematics of Computation*, 31:629-651, 1977.

[13] J. L. García Zapata and J. C. Díaz Martín. A geometric algorithm for winding number computation with complexity analysis. *Journal of Complexity*, 28(3):320-345, 2012.

[14] J. L. García Zapata and J. C. Díaz Martín. Finding the number of roots of a polynomial in a plane region using the winding number. *Computers & Mathematics with Applications*, 67(3):555-568, 2014.

[15] C. Gasquet and P. Witomski. *Fourier Analysis and Applications: Filtering, Numerical Computation, Wavelets*, volume 30. Springer Science & Business Media, 2013.

[16] M. Goldberg. On a boundary extrapolation theorem by Kreiss. *Mathematics of Computation*, 31(138):469-477, 1977.

[17] B. Gustafsson. *High Order Difference Methods for Time Dependent PDE*. Number 38 in Springer series in computational mathematics. Springer, Berlin, 2008.

[18] B. Gustafsson, H.-O. Kreiss, and J. Oliger. *Time-dependent problems and difference methods*. Pure and Applied Mathematics: A Wiley Series of Texts, Monographs and Tracts. Wiley, 2013.

[19] B. Gustafsson, H.-O. Kreiss, and A. Sundström. Stability theory of difference approximations for mixed initial boundary value problems. II. *Mathematics of Computation*, 26(119):649-649, 1972.

[20] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant. Array programming with NumPy. *Nature*, 585(7825):357-362, Sept. 2020.

[21] C. Hermite. Sur la formule d'interpolation de Lagrange. (Extrait d'une lettre de M. Ch. Hermite à M. Borchardt). *Journal für die reine und angewandte Mathematik*, 84:70-79, 1877.

[22] R. Hersh. Mixed problems in several variables. *Journal of Mathematics and Mechanics*, 12(3):317-334, 1963.

[23] H.-O. Kreiss. Stability theory for difference approximations of mixed initial boundary value problems. i. *Mathematics of Computation*, 22(104):703-703, 1968.

[24] S. Lang. *Complex analysis*, volume 103. Springer Science & Business Media, 2013.

[25] P. D. Lax and R. D. Richtmyer. Survey of the stability of linear finite difference equations. *Communications on Pure and Applied Mathematics*, 9(2):267-293, 1956.

[26] T. Li, J. Lu, and C.-W. Shu. Stability analysis of inverse Lax-Wendroff boundary treatment of high order compact difference schemes for parabolic equations. *Journal of Computational and Applied Mathematics*, 400:113711, 2022.

[27] T. Li, C.-W. Shu, and M. Zhang. Stability analysis of the inverse Lax-Wendroff boundary treatment for high order upwind-biased finite difference schemes. *Journal of Computational and Applied Mathematics*, 299:140-158, 2016.

[28] T. Li, C.-W. Shu, and M. Zhang. Stability analysis of the inverse Lax-Wendroff boundary treatment for high order central difference schemes for diffusion equations. *Journal of Scientific Computing*, 70(2):576-607, 2017.

[29] F. Lörcher and C.-D. Munz. Lax-Wendroff-type schemes of arbitrary order in several space dimensions. *IMA Journal of Numerical Analysis*, 27(3):593-615, 2006.

[30] G. Métivier and K. Zumbrun. Symmetrizers and continuity of stable subspaces for parabolic-hyperbolic boundary value problems. *Discrete & Continuous Dynamical Systems - A*, 11(1):205-220, 2004.

[31] A. Meurer, C. P. Smith, M. Paprocki, O. Čertík, S. B. Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J. K. Moore, S. Singh, T. Rathnayake, S. Vig, B. E. Granger, R. P. Muller, F. Bonazzi, H. Gupta, S. Vats, F. Johansson, F. Pedregosa, M. J. Curry, A. R. Terrel, v. Roučka, A. Saboo, I. Fernando, S. Kulal, R. Cimrman, and A. Scopatz. Sympy: symbolic computing in Python. *PeerJ Computer Science*, 3:e103, Jan. 2017.

[32] M. Thuné. Automatic GKS stability analysis. *SIAM J. Sci. Statist. Comput.*, 7(3):959-977, 1986.

[33] F. Vilar and C.-W. Shu. Development and stability analysis of the inverse Lax-Wendroff boundary treatment for central compact schemes. *ESAIM: Mathematical Modelling and Numerical Analysis*, 49(1):39-67, 2015.

[34] L. Wu. The semigroup stability of the difference approximations for initial-boundary value problems. *Mathematics of Computation*, 64(209):71-71, 1995.

UNIV RENNES, CNRS, IRMAR - UMR 6625, F-35000 RENNES, FRANCE.

Email address: benjamin.boutin@univ-rennes1.fr

Email address: pierre.lebarbenchon@univ-rennes1.fr

IMAG, INRIA D'UNIVERSITÉ CÔTE D'AZUR, UNIV. MONTPELLIER, CNRS, MONTPELLIER, FRANCE

Email address: nicolas.seguin@inria.fr
By “constructible function”, we mean here that we establish a computable algorithm to get the matrix $\tilde{\mathcal{B}}(\sigma_0(z), \dots, \sigma_{r-1}(z))$. This algorithm, based on a Gaussian elimination, is fully described in the following proof.

*Proof.* — For $z \in \bar{\mathcal{U}}$ and $\zeta_0 = \sigma_0(z), \dots, \sigma_{r-1} = \sigma_{r-1}(z)$. By a descending induction on $j$ between $m-1$ to $r-1$, we construct a matrix $\mathcal{B}_j(\zeta_0, \dots, \zeta_{r-1}) \in \mathcal{M}_{r,j+1}(\mathbb{C})$ such that

$$ \mathcal{B}\tilde{\mathcal{U}}_{m-1}(z) = \mathcal{B}_j(\zeta_0, \dots, \zeta_{r-1})\tilde{\mathcal{U}}_j(z). $$

*Initialization.* — If $j = m-1$ then one can take $\mathcal{B} = \mathcal{B}_{m-1}(\zeta_0, \dots, \zeta_{r-1})$.

*Induction.* — We assume the induction hypotheses for some $j \in [m-1 : r]$ and we want to prove the result for $j-1$. By equation (4.1), we have $\tilde{\mathcal{U}}_j(z) = \mathcal{P}_j\tilde{\mathcal{U}}_{j-1}(z)$ where

$$ \mathcal{P}_j = \begin{pmatrix} 1 & & & \\ & \ddots & & \\ & & \ddots & \\ & & & 1 \\ (0) & -\zeta_0 & \cdots & -\zeta_{r-1} \end{pmatrix} \in \mathcal{M}_{j+1,j}(\mathbb{C}). $$

We define $\mathcal{B}_{j-1}(\zeta_0, \dots, \zeta_{r-1}) = \mathcal{B}_j(\zeta_0, \dots, \zeta_{r-1})\mathcal{P}_j \in \mathcal{M}_{r,j}(\mathbb{C})$ then we have

$$ \mathcal{B}_{j-1}\tilde{\mathcal{U}}_{j-1}(z) = \mathcal{B}_j\mathcal{P}_j\tilde{\mathcal{U}}_{j-1}(z) = \mathcal{B}_j\tilde{\mathcal{U}}_j(z) = \mathcal{B}\tilde{\mathcal{U}}_{m-1}(z). $$

*Conclusion.* — We define $\tilde{\mathcal{B}}$ by $\mathcal{B}_{r-1}$.

The function $\tilde{\mathcal{B}}$ is easily computable because $(\mathcal{P}_j)_j$ are just matrices of Gaussian elimination. $\square$

By (4.2), the intrinsic Kreiss–Lopatinskii determinant can be written

$$ \begin{aligned} \Delta(z) &= \frac{\det(zK_{0,r-1}(z) - \mathcal{B}K_{0,m-1}(z))}{\det K_{0,r-1}(z)} \\ &= \frac{\det(zK_{0,r-1}(z) - \tilde{\mathcal{B}}K_{0,r-1}(z))}{\det K_{0,r-1}(z)} = \det(zI_r - \tilde{\mathcal{B}}). \quad (4.3) \end{aligned} $$

The matrix $\tilde{\mathcal{B}}$ from Proposition 4.2 depends on coefficients $(\sigma_j(z))_j$. By (4.2), we have

$$ \tilde{\mathcal{B}}(\sigma_0(z), \dots, \sigma_{r-1}(z)) = \mathcal{B}K_{0,m-1}(z)K_{0,r-1}^{-1}(z). $$

Using any computer algebra system, we can compute $\tilde{\mathcal{B}}(\zeta_0, \dots, \zeta_{r-1})$ from the matrix $\mathcal{B}$ and then we compute the coefficients $(\sigma_j(z))_j$ and replace $\zeta_j$ by $\sigma_j(z)$ for all $j \in [0 : r-1]$. It provides that the computation of $\tilde{\mathcal{B}}(\zeta_0, \dots, \zeta_{r-1})$ can be done only once and then applied for different $z$ (and so different $(\sigma_j(z))_j$).
With (4.3), we find again the holomorphic property of $\Delta$ by Lemma 3.1 which states that $z \mapsto \sigma_j(z)$ is holomorphic on $\mathcal{U}$ for all $j \in [0 : r - 1]$.

## 4.2. Computation of $\Delta(\mathcal{S})$

Let us fix a $z_0 \in \mathcal{S}$. To compute $(\sigma_j(z_0))_j$, we need the $r$ roots $(\kappa_j(z_0))_j$ that come from the inside of the unit disk, see Remark 2.7. By the continuity of the roots of polynomial $P_{z_0}$ defined in (3.1) with respect to the parameter $z_0$, for each $\kappa_0(z_0)$ of multiplicity $\beta$ on the unit circle, for a sufficiently small $\varepsilon > 0$, there exists $\eta > 0$ such that for all $z \in B(z_0, \eta)$, the polynomial $P_z$ has exactly $\beta$ roots with multiplicity in $B(\kappa_0(z_0), \varepsilon)$. The explicit value of $\eta$ is given in the following statement.

**LEMMA 4.3.** — Let $z_0$ be on the unit circle. Let $\kappa_0(z_0) \in \mathcal{S}$ be a root of multiplicity $\beta$ of the polynomial $P_{z_0}$ defined in (3.1). Let $\varepsilon > 0$ be such that $\kappa_0(z_0)$ is the only root of $P_{z_0}$ in $\overline{B(\kappa_0(z_0), \varepsilon)}$ and set

$$ \eta = (1 + \varepsilon)^{-r} \min_{\kappa \in \partial B(\kappa_0(z_0), \varepsilon)} |P_{z_0}(\kappa)|. $$

Then for all $z \in B(z_0, \eta)$, the polynomial $P_z$ has exactly $\beta$ roots with multiplicity in $B(\kappa_0(z_0), \varepsilon)$.

The proof of Lemma 4.3 is a consequence of Rouché's theorem, comparing the number of zeros between $P_{z_0}$ and $P_z$ for $z$ close to $z_0$ in $B(\kappa_0(z_0), \varepsilon)$, the details of the proof are not given here, but we refer to [29] for a proof of Rouché's theorem. This lemma is illustrated in Figure 4.1 where $\kappa_0(z_0)$ is of multiplicity 3. The black points are related to $z_0$ and $(\kappa_j(z_0))_j$, and the gray point are related to $z$ and $(\kappa_j(z))_j$.

From the numerical point of view, for a multiple root $\kappa_0(z_0)$, one can take the smallest distance between two roots of $P_{z_0}$ as $\varepsilon$, take $z = (1 + \frac{\eta}{2})z_0 \in \mathcal{U}$. The value of $\eta$ is obtained discretizing the circle of radius $\varepsilon$ centered in $\kappa_0(z_0)$. By Lemma 2.1 (Hersh), there is no roots of $P_z$ on the unit circle, then one can count the roots in $B(\kappa_0(z_0), \varepsilon) \cap \mathbb{D}$ and $B(\kappa_0(z_0), \varepsilon) \cap \mathcal{U}$ to know the number of roots linked to $\kappa_0(z)$ that come from the inside and the outside of the unit disk. After selecting the roots $(\kappa_j(z))_j$ that come from the inside of the unit disk, one may compute their symmetric functions $(\sigma_j(z))_j$. By replacing the formal variables $(\zeta_j)_j$ of $\tilde{B}(\zeta_0, \dots, \zeta_{r-1})$ with $(\sigma_j(z))_j$, one may compute $\Delta(z)$ with the expression (4.3). Instead of computing $(\kappa_j(z))_j$ for each $z$ on the unit circle independently, one may use the continuity of $(\kappa_j(z))_j$ with respect to $z$ in order to describe the movement of the roots $(\kappa_j(z))_j$ for $z \in \mathcal{S}$. After drawing the Kreiss-Lopatinskii curve, the winding number has to be computed in order to use Method 2.14. To do so, we use
<figure><img src="image_8.png" /><figcaption>Figure 4.1. Illustration of Lemma 4.3</figcaption></figure>

the geometric algorithm proposed by García Zapata and Díaz Martín in [15] and [16].

### 4.3. Boundary condition: reconstruction procedure

To define the boundary condition, we use the reconstruction procedure explained in [11]. The framework is the advection equation with a misalignment between the space boundary and the discrete grid points

$$ egin{cases} \partial_t u + a \partial_x u = 0, & t \ge 0, x \in [x_\sigma, 1], \\ u(t, x_\sigma) = g(t), & t \ge 0, \\ u(0, x) = f(x), & x \in [x_\sigma, 1]. \end{cases} \quad (4.4) $$

Without loss of generality, we can assume that $x_\sigma = \sigma \Delta x$ with $\sigma \in [-\frac{1}{2}, \frac{1}{2}]$ (as it is explained in [5]). Let us introduce $x_j$ for $j \Delta x$ and $t^n$ for $n \Delta t$ when $j \ge -r$ and $n \ge 0$. Let $n \in \mathbb{N}$ be a fixed time. The solution $u$ of (4.4) (assumed here to be smooth enough) satisfies

$$ \begin{aligned} & \frac{1}{\Delta x} \int_{x_j - \frac{\Delta x}{2}}^{x_j + \frac{\Delta x}{2}} u(t^n, y) dy \\ & = \frac{1}{\Delta x} \int_{x_j - \frac{\Delta x}{2}}^{x_j + \frac{\Delta x}{2}} \sum_{k=0}^{d-1} \partial_x^k u(t^n, x_\sigma) \frac{(y - x_\sigma)^k}{k!} dy + O(\Delta x^d) \quad (4.5) \end{aligned} $$

using a Taylor expansion of order $d$. Now, let us take a solution $(U_j^n)_{j \ge 0}$ of a scheme of the form (1.2) approximating $u$. Using (4.5), we want to define
the r ghost points $(U_j^n)_{-r \le j \le -1}$. The approximation of equation (4.5) reads, for $j \ge -r$,

$$U_j^n \approx \sum_{k=0}^{d-1} \partial_x^k u(t^n, x_\sigma) \left( \frac{(j + \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} - \frac{(j - \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} \right). \quad (4.6)$$

On the one hand, we use the PDE to convert space derivatives into time derivatives until an index $k_d < d$. The index $k_d$ allows us to know only the first derivatives of the boundary datum $g$ and use extrapolation for the rest. For the advection equation, we have, for all $k \le k_d$,

$$\partial_x^k u(t^n, x_\sigma) = (-a)^{-k} \partial_t^k u(t^n, x_\sigma) = (-a)^{-k} g^{(k)}(t^n).$$

Equation (4.6) becomes, for $j \ge -r$,

$$U_j^n \approx \sum_{k=0}^{k_d} (-a)^{-k} g^{(k)}(t^n) \left( \frac{(j + \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} - \frac{(j - \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} \right) + \sum_{k=k_d+1}^{d-1} \partial_x^k u(t^n, x_\sigma) \left( \frac{(j + \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} - \frac{(j - \frac{1}{2} - \sigma)^{k+1}}{(k+1)!} \right). \quad (4.7)$$

On the other hand, we need to define $(\partial_x^k u(t^n, x_\sigma))_{k=k_d+1}^{d-1}$, but using (4.7) for $j \in [0 : d - k_d - 2]$, we can deduce the unknowns $(U_j^n)_{-r \le j \le -1}$. Writing

$$\mathfrak{A}_- = (U_{-r}^n, \dots, U_{-1}^n)^T, \quad \mathfrak{A}_+ = (U_0^n, \dots, U_{d-k_d-2}^n)^T \\ \text{and} \quad \Theta^n = (\partial_x^{k_d+1} u(t^n, x_\sigma), \dots, \partial_x^{d-1} u(t^n, x_\sigma))^T,$$

we have a condensed formulation of (4.7):

$$\left\{ \begin{array}{l} \mathfrak{A}_- = S_-^n + \mathcal{Y}_- \Theta^n, \\ \mathfrak{A}_+ = S_+^n + \mathcal{Y}_+ \Theta^n, \end{array} \right. \qquad (4.8)$$

where $S_-^n \in \mathbb{R}^r$, $S_+^n \in \mathbb{R}^{d-k_d-1}$, $\mathcal{Y}_- \in \mathcal{M}_{r,d-k_d-1}(\mathbb{R})$ and $\mathcal{Y}_+ \in \mathcal{M}_{d-k_d-1}(\mathbb{R})$ with

$$(S_-^n)_i = \sum_{k=0}^{k_d} (-a)^{-k} g^{(k)}(t^n) \left( \frac{\left(-i + \frac{1}{2} - \sigma\right)^{k+1}}{(k+1)!} - \frac{\left(-i - \frac{1}{2} - \sigma\right)^{k+1}}{(k+1)!} \right), \\ (S_+^n)_i = \sum_{k=0}^{k_d} (-a)^{-k} g^{(k)}(t^n) \left( \frac{\left(i - \frac{1}{2} - \sigma\right)^{k+1}}{(k+1)!} - \frac{\left(i - 1 - \frac{1}{2} - \sigma\right)^{k+1}}{(k+1)!} \right),$$
$$ (\mathcal{Y}_-)_i, j = \left( \frac{\left( -(r-i+1) + \frac{1}{2} - \sigma \right)^{j+k_d+1}}{(j+k_d+1)!} - \frac{\left( -(r-i+1) - \frac{1}{2} - \sigma \right)^{j+k_d+1}}{(j+k_d+1)!} \right). $$

$$ (\mathcal{Y}_+)_i, j = \left( \frac{\left( i-1+\frac{1}{2}-\sigma \right)^{j+k_d+1}}{(j+k_d+1)!} - \frac{\left( i-1-\frac{1}{2}-\sigma \right)^{j+k_d+1}}{(j+k_d+1)!} \right), $$

Eliminating the space derivatives $\Theta^n$ in (4.8) gives the following boundary condition:

$$ \mathfrak{U}_- = \mathcal{Y}_- \mathcal{Y}_+^{-1} \mathfrak{U}_+ + \mathcal{S}_-^n - \mathcal{Y}_- \mathcal{Y}_+^{-1} \mathcal{S}_+^n. \quad (4.9) $$

Equation (4.9) is exactly the boundary equations (1.5) of the scheme which define the r ghost points of the scheme. To write the boundary condition as equation (1.7) with expression (1.8), we identify

$$ B \stackrel{\text{def}}{=} \mathcal{Y}_- \mathcal{Y}_+^{-1} \text{ and } \begin{pmatrix} g_-^n \\ \vdots \\ g_-^n \end{pmatrix} \stackrel{\text{def}}{=} \mathcal{S}_-^n - \mathcal{Y}_- \mathcal{Y}_+^{-1} \mathcal{S}_+^n. $$

As in [11], $\mathcal{R}^{d,k_d}$ denotes the reconstruction procedure where $d$ is the order of consistency of the method and $k_d$ the index when we change from time derivatives to extrapolation. For example, the reconstruction procedure $\mathcal{R}^{3,0}$ for $r=2$ and $\sigma=0.4$ leads to

$$ \mathcal{Y}_- = \begin{pmatrix} -\frac{12}{5} & \frac{1753}{600} \\ -\frac{7}{5} & \frac{613}{600} \end{pmatrix}, \quad \mathcal{Y}_+ = \begin{pmatrix} -\frac{2}{5} & \frac{73}{600} \\ \frac{3}{5} & \frac{133}{600} \end{pmatrix} \quad \text{and} \quad B = \begin{pmatrix} \frac{1371}{97} & \frac{526}{97} \\ \frac{554}{97} & \frac{143}{97} \end{pmatrix} $$

### 4.4. Example of O3 scheme

As it is done in [11], we want to find the stability area for the O3 scheme defined, for $j \in \mathbb{N}$ and $n \in \mathbb{N}$, by

$$ U_j^{n+1} = \left( \frac{\lambda^3}{6} - \frac{\lambda}{6} \right) U_{j-2}^n + \left( \lambda + \frac{\lambda^2}{2} - \frac{\lambda^3}{2} \right) U_{j-1}^n \\ + \left( 1 - \frac{\lambda}{2} - \lambda^2 + \frac{\lambda^3}{2} \right) U_j^n + \left( \frac{\lambda^2}{2} - \frac{\lambda^3}{6} - \frac{\lambda}{3} \right) U_{j+1}^n \quad (4.10) $$

The O3 scheme is a scheme with $r=2$ and $p=1$ and is Cauchy-stable for $\lambda \in ]0, 1]$ and is of order 3 both in time and space. The reconstruction $\mathcal{R}^{3,0}$ for the O3 scheme and $\sigma=0.4$ leads to

$$ B = \begin{pmatrix} \frac{1371}{97} & \frac{526}{97} & 0 \\ \frac{554}{97} & \frac{143}{97} & 0 \end{pmatrix} \text{ and} $$

$$ B = \begin{pmatrix} \frac{180\lambda^2}{97} + \frac{277\lambda}{97} + 1 & \frac{120\lambda^2}{97} + \frac{23\lambda}{97} & 0 \\ \frac{263\lambda^3}{582} + \frac{\lambda^2}{2} + \frac{14\lambda}{291} & \frac{217\lambda^3}{291} - \lambda^2 - \frac{217\lambda}{291} + 1 & -\frac{\lambda^3}{6} + \frac{\lambda^2}{2} - \frac{\lambda}{3} \end{pmatrix}. $$
Using the reformulation (4.3) of the Kreiss–Lopatinskii determinant, we have

$$ \tilde{B}(\zeta_0, \zeta_1) = B \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ -\zeta_0 & -\zeta_1 \end{pmatrix} = \left( \begin{array}{c} \frac{180\lambda^2}{97} + \frac{277\lambda}{97} + 1 \\ \frac{(263+97\zeta_0)\lambda^3}{582} + \frac{(1-\zeta_0)\lambda^2}{2} + \frac{(14+97\zeta_0)\lambda}{291} \\ \frac{434+97\zeta_1)\lambda^3}{582} - \frac{120\lambda^2}{97} + \frac{23\lambda}{97} - \frac{(2+\zeta_1)\lambda^2}{2} - \frac{(217-97\zeta_1)\lambda}{291} + 1 \end{array} \right). $$

For example, for $\sigma = 0.4$, Figure 4.2 shows that the O3 scheme with $\mathcal{R}^{3,0}$ boundary is stable for $\lambda = 0.4$ (because $r - \text{Ind}_{\Delta(\mathcal{S})}(0) = 0$) and is unstable for $\lambda = 0.9$ (because $r - \text{Ind}_{\Delta(\mathcal{S})}(0) = 1$).

<figure><img src="image_9.png" /><figcaption>Figure 4.2. Curve $\Delta(\mathcal{S})$ for O3 scheme for $\sigma = 0.4$, for $\lambda \in \{0.4, 0.9\}$ with reconstruction boundary $\mathcal{R}^{3,0}$.</figcaption></figure>

We can draw the same figure as the Figure 4 of [11] but instead of using a computation of the spectral radius of the truncated quasi-Toeplitz matrix, we use our strategy of counting the number of instability modes, see Figure 4.3 which is much more reliable (since it is parameter free) and efficient. In Figure 4.3, every area stamped with 0 is a domain where the O3 scheme is stable. The odd pattern for very small $\lambda$ (approximatively between 0 and 0.01) of Figure 4.3 may be due to difficulties for computing the winding number. Indeed, for very small values of $\lambda$, the Kreiss–Lopatinskii determinant is really close to the origin and even with a refinement (see the next example for more details on this procedure), the computation of the winding number may become inaccurate, which is not a problem in practice since it would correspond to very small, then unusable, time steps.
<figure><img src="image_10.png" /><figcaption>Figure 4.3. Number of zeros of the Kreiss–Lopatinskii determinant of O3 scheme with different reconstruction boundaries for $\lambda \in ]0, 1]$ and $\sigma \in ]-0.5, 0.5[$.</figcaption></figure>

### 4.5. Example of Lax–Wendroff 5

The fifth-order Lax–Wendroff scheme, called LW5, which has been proposed in [36], can be written, for all $n \in \mathbb{N}$, for all $j \in \mathbb{N}$, as

$$ 
\begin{aligned}
U_{j}^{n+1} = & \frac{\lambda(\lambda-2)(\lambda-1)(\lambda+1)(\lambda+2)}{120} U_{j-3}^{n} - \frac{\lambda(\lambda-1)(\lambda-3)(\lambda+1)(\lambda+2)}{24} U_{j-2}^{n} \\
& + \frac{\lambda(\lambda-2)(\lambda-3)(\lambda+1)(\lambda+2)}{12} U_{j-1}^{n} + \left(1 - \frac{\lambda(\lambda^4-3\lambda^3-5\lambda^2+15\lambda+4)}{12}\right) U_{j}^{n} \\
& + \frac{\lambda(\lambda-1)(\lambda-2)(\lambda-3)(\lambda+2)}{24} U_{j+1}^{n} - \frac{\lambda(\lambda-1)(\lambda-2)(\lambda-3)(\lambda+1)}{120} U_{j+2}^{n}.
\end{aligned} 
\quad (4.11) $$

This scheme LW5 is Cauchy-stable for $\lambda \in ]0, 1]$ and is of order 5 both in time and space.

Figure 4.4 illustrates the computation of the number of instabilities for LW5 for different reconstruction boundaries where $\sigma = 0.4$ with respect to $\lambda \in ]0, 1]$. As in the previous example, it may happen that the Kreiss–Lopatinskii curve is too close to the origin, the winding number of the origin cannot be correctly computed. Following the geometric algorithm proposed by García Zapata in [16], a refinement of the discretization then improves the effective computation of the winding number. Figure 4.5 represents such refinement with close-up close to the origin. However, even with this strategy,
<figure><img src="image_11.png" /><figcaption>Figure 4.4. Number of zeros of the Kreiss–Lopatinskii determinant of LW5 scheme with different reconstruction boundaries for $\lambda \in ]0, 1]$ and $\sigma = 0.4$.</figcaption></figure>

for very small values of $\lambda$, we cannot refine more than the machine precision, that is why there is still some odd pattern for very small $\lambda$ in Figure 4.3, Figure 4.4 and Figure 4.6. As we already discussed, such very small time step are however not used in practice. The stability area with respect to both parameters $\lambda$ and $\sigma$ are drawn in Figure 4.6, considering again successively various reconstruction boundary conditions.

All the figures can be easily computed in Python with the common NumPy [24] library and the SymPy [38] library for the computer algebra system. The algorithm is really efficient. For each subfigure of Figure 4.3, the 1600 runs takes less than a couple of minutes of computation achieved on a standard laptop. Moreover, our procedure provides sharp results. In particular, contrary to numerical investigations of stability which are based on the computation of the spectral radius, no arbitrary truncation of (quasi-)Toeplitz matrices is needed.

## 5. Future directions

The present theoretical and numerical results are restricted to the class of one-step schemes. Actually, the multistep case presents similarities with the one-time step case thanks to the work of Coulombel [7]. Hence, Theorem 2.11 still holds and Corollary 2.13 has to be adapted: the number of zeros in $\mathcal{U}$ of
<figure><img src="image_12.png" /><figcaption>Figure 4.5. (a) representation of Δ(S) for LW5 with the boundary condition $\mathcal{R}^{6,1}$, for $\lambda = 0.01$ and $\sigma = 0$, zoom (b) without refinement (c) and with refinement (d).</figcaption></figure>

the equation $\Delta(z) = 0$ is now $r(s+1) - \text{Ind}_{\Delta(\mathcal{S})}(0)$ where $s+1$ is the number of time-steps of the scheme, this work is initiated in [32, Sec. 6.7]. In another direction, for implicit schemes or for more general boundary conditions, such as absorbing boundary conditions [13] and [14] or transparent boundary conditions [1] and [8], it seems to be more challenging to have a such easy-to-use theory. The work of Benoit [3] may be a source of inspiration to work on a bounded space domain. In that case, we have to modify the stability estimates (1.9) to deal with the two boundaries and introduce two Kreiss-Lopatinskii determinants, one for each side.
<figure><img src="image_13.png" /><figcaption>Figure 4.6. Number of zeros of the Kreiss–Lopatinskii determinant of LW5 scheme with different reconstruction boundaries for $\lambda \in ]0, 1]$ and $\sigma \in ]-0.5, 0.5[$.</figcaption></figure>
Bibliography

[1] A. ARNOLD, M. EHRHARDT & I. SOFRONOV, "Discrete transparent boundary conditions for the Schrödinger equation: fast calculation, approximation, and stability", Commun. Math. Sci. 1 (2003), no. 3, p. 501-556.

[2] R. M. BEAM & R. F. WARMING, "The asymptotic spectra of banded Toeplitz and quasi-Toeplitz matrices", *SIAM J. Sci. Comput.* **14** (1993), no. 4, p. 971-1006.

[3] A. BENOIT, "Stability of finite difference schemes approximation for hyperbolic boundary value problems in an interval", *Math. Comput.* **91** (2022), no. 335, p. 1171-1212.

[4] S. BENZONI-GAVAGE & D. SERRE, *Multidimensional hyperbolic partial differential equations: first-order systems and applications*, Oxford Mathematical Monographs, Clarendon Press, 2007.

[5] B. BOUTIN, P. LE BARBENCHON & N. SEGUIN, "On the stability of totally upwind schemes for the hyperbolic initial boundary value problem", *IMA J. Numer. Anal.* **44** (2024), no. 2, p. 1211-1241.

[6] B. BOUTIN, T. H. T. NGUYEN, A. SYLLA, S. TRAN-TIEN & J.-F. COULOMBEL, "High order numerical schemes for transport equations on bounded domains", *ESAIM: ProcS* **70** (2021), p. 84-106.

[7] J.-F. COULOMBEL, "Stability of finite difference schemes for hyperbolic initial boundary value problems", in *HCDTE lecture notes. Part I. Nonlinear hyperbolic PDEs, dispersive and transport equations*, AIMS Series on Applied Mathematics, vol. 6, American Institute of Mathematical Sciences, 2013, p. 98-225.

[8] ——, "Transparent numerical boundary conditions for evolution equations: Derivation and stability analysis", *Ann. Fac. Sci. Toulouse, Math.* **28** (2019), no. 2, p. 259-327.

[9] R. COURANT, K. FRIEDRICHS & H. LEWY, "Über die partiellen differenzengleichungen der mathematischen physik", *Math. Ann.* **100** (1928), no. 1, p. 32-74.

[10] J. CRANK & P. NICOLSON, "A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type", *Proc. Camb. Philos. Soc.* **43** (1947), p. 50-67.

[11] G. DAKIN, B. DESPRÉS & S. JAOUEN, "Inverse Lax-Wendroff boundary treatment for compressible lagrange-remap hydrodynamics on cartesian grids", *J. Comput. Phys.* **353** (2018), p. 228-257.

[12] S. DEL PINO & H. JOURDREN, "Arbitrary high-order schemes for the linear advection and wave equations: application to hydrodynamics and aeroacoustics", *C. R. Math.* **342** (2006), no. 6, p. 441-446.

[13] M. EHRHARDT, "Absorbing boundary conditions for hyperbolic systems", *Numer. Math., Theory Methods Appl.* **3** (2010), no. 3, p. 295-337.

[14] B. ENGQUIST & A. MAJDA, "Absorbing boundary conditions for the numerical simulation of waves", *Math. Comput.* **31** (1977), no. 139, p. 629-651.

[15] J.-L. GARCÍA ZAPATA & J. C. DÍAZ MARTÍN, "A geometric algorithm for winding number computation with complexity analysis", *J. Complexity* **28** (2012), no. 3, p. 320-345.

[16] ——, "Finding the number of roots of a polynomial in a plane region using the winding number", *Comput. Math. Appl.* **67** (2014), no. 3, p. 555-568.

[17] C. GASQUET & P. WITOMSKI, *Fourier analysis and applications. Filtering, numerical computation, wavelets. Transl. from the French by R. Ryan*, Texts in Applied Mathematics, vol. 30, Springer, 1999.

[18] M. GOLDBERG, "On a boundary extrapolation theorem by Kreiss", *Math. Comput.* **31** (1977), no. 138, p. 469-477.
[19] M. GOLDBERG & E. TADMOR, "Scheme-independent stability criteria for difference approximations of hyperbolic initial-boundary value problems. I", *Math. Comput.* **32** (1978), no. 144, p. 1097-1107.

[20] ——, "Scheme-independent stability criteria for difference approximations of hyperbolic initial-boundary value problems. II", *Math. Comput.* **36** (1981), no. 154, p. 603-626.

[21] B. GUSTAFSSON, *High Order Difference Methods for Time Dependent PDE*, Springer Series in Computational Mathematics, vol. 38, Springer, 2008.

[22] B. GUSTAFSSON, H.-O. KREISS & J. OLIGER, *Time-dependent problems and difference methods*, 2nd ed., Pure and Applied Mathematics, John Wiley & Sons, 2013.

[23] B. GUSTAFSSON, H.-O. KREISS & A. SUNDSTRÖM, "Stability theory of difference approximations for mixed initial boundary value problems. II", *Math. Comput.* **26** (1972), no. 119, p. 649-686.

[24] C. R. HARRIS, K. J. MILLMAN, S. J. VAN DER WALT, R. GOMMERS, P. VIRTANEN, D. COURNAPEAU, E. WIESER, J. TAYLOR, S. BERG, N. J. SMITH, R. KERN, M. PICUS, S. HOYER, M. H. VAN KERKIJK, M. BRETT, A. HALDANE, J. F. DEL RÍO, M. WIEBE, P. PETERSON, P. GÉRARD-MARCHANT, K. SHEPPARD, T. REDDY, W. WECKESSER, H. ABBASI, C. GOHLKE & T. E. OLIPHANT, "Array programming with NumPy", *Nature* **585** (2020), no. 7825, p. 357-362.

[25] M. C. HERMITE & M. BORCHARDT, "Sur la formule d'interpolation de Lagrange", *J. Reine Angew. Math.* **84** (1878), p. 70-79.

[26] R. HERSH, "Mixed problems in several variables", *J. Math. Mech.* **12** (1963), p. 317-334.

[27] A. ISERLES & G. STRANG, "The optimal accuracy of difference schemes", *Trans. Am. Math. Soc.* **277** (1983), no. 2, p. 779-803.

[28] H.-O. KREISS, "Stability theory for difference approximations of mixed initial boundary value problems. I", *Math. Comput.* **22** (1968), p. 703-714.

[29] S. LANG, *Complex analysis*, 4 ed., Graduate Texts in Mathematics, vol. 103, Springer, 1999.

[30] P. D. LAX & R. D. RICHTMYER, "Survey of the stability of linear finite difference equations", *Commun. Pure Appl. Math.* **9** (1956), no. 2, p. 267-293.

[31] P. LE BARBENCHON, "Boundaryscheme: package Python for numerical schemes with boundaries", 2023, https://zenodo.org/doi/10.5281/zenodo.7773741.

[32] ——, "Étude théorique et numérique de la stabilité GKS pour des schémas d'ordre élevé en présence de bords", PhD Thesis, Université de Rennes, 2023, https://theses.hal.science/tel-04214887.

[33] T. LI, J. LU & C.-W. SHU, "Stability analysis of inverse Lax-Wendroff boundary treatment of high order compact difference schemes for parabolic equations", *J. Comput. Appl. Math.* **400** (2022), article no. 113711 (26 pages).

[34] T. LI, C.-W. SHU & M. ZHANG, "Stability analysis of the inverse Lax-Wendroff boundary treatment for high order upwind-biased finite difference schemes", *J. Comput. Appl. Math.* **299** (2016), p. 140-158.

[35] ——, "Stability analysis of the inverse Lax-Wendroff boundary treatment for high order central difference schemes for diffusion equations", *J. Sci. Comput.* **70** (2017), no. 2, p. 576-607.

[36] F. LÖRCHER & C.-D. MUNZ, "Lax-Wendroff-type schemes of arbitrary order in several space dimensions", *IMA J. Numer. Anal.* **27** (2006), no. 3, p. 593-615.

[37] G. MÉTIVIER & K. ZUMBRUN, "Symmetrizers and continuity of stable subspaces for parabolic-hyperbolic boundary value problems", *Discrete Contin. Dyn. Syst.* **11** (2004), no. 1, p. 205-220.
[38] A. MEURER, C. P. SMITH, M. PAPROCKI, O. ČERTÍK, S. B. KIRPICHEV, M. ROCKLIN, A. KUMAR, S. IVANOV, J. K. MOORE, S. SINGH, T. RATHNAYAKE, S. VIG, B. E. GRANGER, R. P. MULLER, F. BONAZZI, H. GUPTA, S. VATS, F. JOHANSSON, F. PEDREGOSA, M. J. CURRY, A. R. TERREL, S. ROUČKA, A. SABOO, I. FERNANDO, S. KULAL, R. CIMRMAN & A. SCOPATZ, "SymPy: symbolic computing in Python", *PeerJ Comput. Sci.* **3** (2017), article no. e103.

[39] M. THUNÉ, "Automatic GKS stability analysis", *SIAM J. Sci. Stat. Comput.* **7** (1986), no. 3, p. 959-977.

[40] F. VILAR & C.-W. SHU, "Development and stability analysis of the inverse Lax-Wendroff boundary treatment for central compact schemes", *ESAIM, Math. Model. Numer. Anal.* **49** (2015), no. 1, p. 39-67.

[41] L. WU, "The semigroup stability of the difference approximations for initial-boundary value problems", *Math. Comput.* **64** (1995), no. 209, p. 71-88.