# Amplification of numerical wave packets for transport equations with two boundaries

Romain Bonnet-Eymard, Jean-François Coulombel, Grégory Faye

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

## To cite this version:

Romain Bonnet-Eymard, Jean-François Coulombel, Grégory Faye. Amplification of numerical wave packets for transport equations with two boundaries. CNRS; ENS Paris Saclay. 2025. hal-05013751

HAL Id: hal-05013751

https://hal.science/hal-05013751v1

Submitted on 31 Mar 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Amplification of numerical wave packets for transport equations with two boundaries

Romain BONNET-EYMARD*, Jean-François COULOMBEL & Grégory FAYE†

March 31, 2025

## Abstract

The purpose of this note is to investigate the coupling of Dirichlet and Neumann numerical boundary conditions for the transport equation set on an interval. When one starts with a stable finite difference scheme on the lattice $Z$ and each numerical boundary condition is taken separately with the Neumann extrapolation condition at the outflow boundary, the corresponding numerical semigroup on a half-line is known to be bounded. It is also known that the coupling of such numerical boundary conditions on a compact interval yields a stable approximation, even though large time exponentially growing modes may occur. We review the different stability estimates associated with these numerical boundary conditions and give explicit examples of such exponential growth phenomena for finite difference schemes with “small” stencils. This provides numerical evidence for the optimality of some stability estimates on the interval.

## 1 Introduction

We consider the following transport problem. Given a positive velocity $a > 0$ and an interval length $L > 0$, we consider the transport equation at velocity $a$ on the interval $(0, L)$ with homogeneous Dirichlet condition at the incoming boundary (that is, at $x = 0$ here since $a$ is positive):

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

It is understood that $f$ is a given real valued function on $(0, L)$ that is, say, at least continuous on the closed interval $[0, L]$ with $f(0) = 0$ so that the compatibility condition at the corner $t = x = 0$ is satisfied. Extending $f$ by 0 on the set $\mathbb{R}^-$ of negative real numbers, the solution to (1) is given by the formula:

$$ \forall t \ge 0, \quad \forall x \in (0, L), \quad u(t, x) = f(x - at). \qquad (2) $$

<footnote>*École Normale Supérieure Paris-Saclay, 4 avenue des Sciences, 91190, Gif-Sur-Yvette, France. Email : rbonnete@ens-paris-saclay.fr</footnote>

<footnote>†Institut de Mathématiques de Toulouse - UMR 5219, Université de Toulouse ; CNRS, Université Paul Sabatier, 118 route de Narbonne, 31062 Toulouse Cedex 9, France. Research of J.-F. C. was supported by ANR project HEAD under grant agreement ANR-24-CE40-3260. G.F. acknowledges support from Labex CIMI under grant agreement ANR-11-LABX-0040, from ANR project Indyana under grant agreement ANR-21-CE40-0008 and from ANR project HEAD under grant agreement ANR-24-CE40-3260. Emails: jean-francois.coulombel@math.univ-toulouse.fr, gregory.faye@math.univ-toulouse.fr</footnote>
In particular, the solution $u(t, \cdot)$ vanishes on $(0, L)$ for $t \ge L/a$.

The goal of this note is to investigate the numerical counterpart of this rather trivial problem. We do not aim at the most general framework and therefore make several simplifying assumptions. First of all, we consider once and for all a fixed parameter $\lambda > 0$. Given the interval length $L > 0$, we consider an integer $J \ge 1$ that is meant to be large, and we define the space step $\Delta x := L/(J + 1)$. The grid points are labeled as $x_j := j \Delta x$ for any $j \in \mathbb{Z}$, in such a way that the interval $(0, L)$ is divided into the cells $(x_j, x_{j+1})$ with $j = 0, \dots, J$. The time step $\Delta t$ is then defined as $\Delta t := \lambda \Delta x$. For a given integer $n \in \mathbb{N}$ and $j \in \{0, \dots, J\}$, we let $u_j^n$ denote the approximation of the solution to (1) within the cell $(x_j, x_{j+1}) \times (n\Delta t, (n+1)\Delta t)$. The considered numerical scheme will update the vector $(u_0^n, \dots, u_J^n)$ into a new vector $(u_0^{n+1}, \dots, u_J^{n+1})$ for each $n \in \mathbb{N}$.

For convenience, we introduce the following notation for the discrete difference operators with respect to the spatial index $j$. Given a collection of three real numbers $u_{j-1}, u_j, u_{j+1}$, the discrete derivatives $(D_-u)_j$ and $(D_+u)_j$ at the index $j$ are defined as:

$$ (D_{-}u)_j := u_j - u_{j-1}, \quad (D_{+}u)_j := u_{j+1} - u_j, $$

and we view $D_-$, $D_+$ as operators acting on either sequences or vectors meaning that we allow ourselves to iterate them. For simplicity, we shall omit most of the time to write the brackets and simply use the notation $D_-u_j$ or $D_+u_j$. Such notation is used, of course, at any integer $j$ for which the action of the operator makes sense. As an example, we have:

$$ D_-^2 u_j = u_j - 2u_{j-1} + u_{j-2}, $$

and larger powers of $D_-$ or $D_+$ can be computed by making appeal to binomial coefficients.

**The considered numerical scheme.** We consider two integers $r, p \in \mathbb{N}$ and some real coefficients $a_{-r}, \dots, a_p$ with $a_{-r} \ne 0$ and $a_p \ne 0$. These coefficients are assumed to depend only on $\lambda$ and $a$. They should satisfy the consistency conditions:

$$ \sum_{\ell=-r}^{p} a_{\ell} = 1, \quad \sum_{\ell=-r}^{p} \ell a_{\ell} = -\lambda a. \qquad (3) $$

The numerical scheme in the interior domain reads:

$$ \forall j = 0, \dots, J, \quad u_j^{n+1} = \sum_{\ell=-r}^{p} a_{\ell} u_{j+\ell}^{n}. \qquad (4) $$

It is understood that $r$ and $p$ are fixed while the number $J$ of cells may get arbitrarily large. The update (4) from the discrete time $n$ to the following time level $n + 1$ is possible only if we have the ghost cell values $u_{-r}^n, \dots, u_{-1}^n$ and $u_{J+1}^n, \dots, u_{J+p}^n$ at our disposal. The ghost cell values on the left of the interval are meant to provide for approximations of the trace of the solution so it seems reasonable, in view of (1), to impose the following homogeneous Dirichlet condition:

$$ \forall n \in \mathbb{N}, \quad \forall \mu = -r, \dots, -1, \quad u_{\mu}^{n} = 0. \qquad (5) $$

This is not the only option but it will be the one we choose here for the sake of simplicity. Prescribing numerical boundary conditions on the right of the interval is a little more tricky since the continuous problem (1) does not impose anything at first glance. We follow here the extrapolation procedure that
has been analyzed in [3, 5, 6] and other works. We thus consider a fixed integer $k \in \mathbb{N}$ and define the ghost cell values $u_{J+1}^n, \dots, u_{J+p}^n$ by imposing:

$$ \forall n \in \mathbb{N}, \quad \forall \mu = 1, \dots, p, \quad D_{-}^{k} u_{J+\mu}^{n} = 0. \qquad (6) $$

The numerical boundary conditions can be understood as follows: the condition (6) determines $u_{J+1}^n$ in terms of interior values $u_{J+1-p}^n, \dots, u_J^n$ that are known in advance (take first $\mu = 1$ in (6)). One then determines iteratively $u_{J+2}^n, \dots, u_{J+p}^n$ as when one solves a lower triangular linear system. As for $r$ and $p$, it is understood that $k$ is fixed while $J$ is large and satisfies at least $J + 1 - k \ge 0$ so that the ghost cell values $u_{J+1}^n, \dots, u_{J+p}^n$ are all determined from (6) by using the interior values $u_0^n, \dots, u_J^n$.

The numerical scheme (4), (5), (6) is initialized by considering:

$$ \forall j = 0, \dots, J, \quad u_j^0 := \frac{1}{\Delta x} \int_{x_j}^{x_{j+1}} f(x) dx, $$

where we recall that $f$ stands for the initial condition in (1) and $x_0 = 0, x_{J+1} = L$. The ghost cell values $u_{-r}^0, \dots, u_{-1}^0$ and $u_{J+1}^0, \dots, u_{J+p}^0$ at the initial time $n = 0$ are determined by (5) and (6).

There are two ways to consider and implement the numerical scheme (4), (5), (6). One can either consider it as a time iteration on the vector $(u_0^n, \dots, u_J^n) \in \mathbb{R}^{J+1}$ and write it as:

$$ \forall n \in \mathbb{N}, \quad \begin{pmatrix} u_0^{n+1} \\ \vdots \\ u_J^{n+1} \end{pmatrix} = A \begin{pmatrix} u_0^n \\ \vdots \\ u_J^n \end{pmatrix}, \qquad (7) $$

for a convenient square matrix $A \in \mathcal{M}_{J+1}(\mathbb{R})$. Or one can consider (4), (5), (6) as a time iteration on the extended vector $(u_{-r}^n, \dots, u_{J+p}^n) \in \mathbb{R}^{J+p+r+1}$ but the set of all possible initial conditions is submitted to the restrictions imposed by (5) or (6) so the stability problem is less easy to rewrite in that framework (though the implementation might be easier). We shall therefore consider here the formulation (7), which amounts to considering all ghost cell values as auxiliary data that are necessary in the iteration process but that are not meaningful in the stability analysis.

The main stability problem that we investigate here is inspired from [1] and amounts to determining whether the iteration matrix $A$ is power bounded:

$$ \sup_{n \in \mathbb{N}} \|A^n\| < +\infty, \qquad (8) $$

where $\| \cdot \|$ is, at this stage, any norm on $\mathcal{M}_{J+1}(\mathbb{R})$ since all norms are equivalent on that space. As is well-known, a necessary condition for (8) to hold is that the spectral radius of $A$ should not be larger than 1. We do not discuss here the choice of the norm and the dependence on the (large) dimension $J$ even though this issue would be extremely meaningful in a discussion about convergence estimates. Our primary focus here is to determine whether the spectral radius of $A$ can exceed 1.

We first review some known stability estimates for the numerical scheme (4), (5), (6) and the two related problems on a half-line. We then present some numerical evidence for the existence of exponentially growing numerical wave packets that may occur even for the very first examples $k = 1$ (first order extrapolation) or $k = 2$ (second order extrapolation). These examples indicate that the stability condition exhibited in [1] is not automatically satisfied in the case of the Dirichlet and Neumann condition when the stencil of the numerical scheme is arbitrary.
# 2 A reminder on Dirichlet and Neumann numerical boundary conditions

## 2.1 Three point schemes

We first recall why the case of three point schemes and $k = 1$ (first order extrapolation) can be easily handled. We consider the case $p = r = 1$, so that the three coefficients $a_{-1}, a_0, a_1$ satisfy (3). The analysis below even allows $p$ (and therefore $a_1$) to be zero. We can equivalently rewrite the iteration (4) as:

$$ \forall j = 0, \dots, J, \quad u_j^{n+1} = u_j^n - \frac{\lambda a}{2} (u_{j+1}^n - u_{j-1}^n) + \frac{\nu}{2} (u_{j+1}^n - 2u_j^n + u_{j-1}^n), \quad (9) $$

where $\nu \in \mathbb{R}$ is a parameter that may depend on $\lambda$ and $a$. Typical choices are:

* $\nu = 1$; one then obtains the so-called Lax-Friedrichs scheme.

* $\nu = \lambda a$; one then obtains the so-called upwind scheme (for which $p = 0$ and $a_1 = 0$).

* $\nu = (\lambda a)^2$; one then obtains the so-called Lax-Wendroff scheme.

The interior scheme (9) is combined with the Dirichlet condition (5) on the left of the interval, that is, $u_{-1}^n = 0$, and the first order extrapolation (or Neumann) condition (6) with $k = 1$ on the right of the interval, that is, $u_{j+1}^n = u_j^n$. Setting:

$$ a_{-1} := \frac{\lambda a + \nu}{2}, \quad a_0 := 1 - \nu, \quad a_1 := \frac{-\lambda a + \nu}{2}, $$

the associated iteration matrix $A$ in (7) reads:

$$ A = \begin{pmatrix} a_0 & a_1 & 0 & \cdots & \cdots & \cdots & 0 \\ a_{-1} & a_0 & a_1 & \cdots & & & & \vdots \\ 0 & \ddots & \ddots & \ddots & \ddots & & & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ & & \ddots & a_{-1} & a_0 & a_1 & 0 & \vdots \\ & & & 0 & a_{-1} & a_0 & a_1 & \vdots \\ 0 & \cdots & \cdots & \cdots & 0 & a_{-1} & a_0 + a_1 & \vdots \end{pmatrix} \in \mathcal{M}_{J+1}(\mathbb{R}). \quad (10) $$

A straightforward result is the following:

**Lemma 1.** Let the parameters in (9) satisfy:

$$ \lambda a \in [0, 1], \quad \nu \in [(\lambda a)^2, 1]. $$

Then the matrix $A$ in (10) is power bounded.

This is, to some extent, the most favorable case where the coupling of Dirichlet and Neumann conditions at the inflow and outflow boundaries yields a straightforward stability estimate.
Proof of Lemma 1. We introduce the standard $\ell^2$-norm on $\mathbb{R}^{J+1}$:

$$ \forall X \in \mathbb{R}^{J+1}, \quad |X|^2 := \sum_{j=0}^{J} X_j^2, $$

and let $\| \cdot \|$ denote the associated matrix norm on $\mathcal{M}_{J+1}(\mathbb{R})$. The argument below will yield $\|A\| \le 1$, showing in particular that $A$ is power bounded since $\| \cdot \|$ is a matrix norm. The bound for the powers of $A$ is even uniform here with respect to the dimension $J$.

We consider a vector $(u_0, \dots, u_J) \in \mathbb{R}^{J+1}$ and define the ghost cell values $u_{-1} := 0$ and $u_{J+1} := u_J$. We then define the vector $(v_0, \dots, v_J) \in \mathbb{R}^{J+1}$ as:

$$ \begin{pmatrix} v_0 \\ \vdots \\ v_J \end{pmatrix} = A \begin{pmatrix} u_0 \\ \vdots \\ u_J \end{pmatrix}, $$

so that we equivalently have:

$$ \forall j = 0, \dots, J, \quad v_j = a_{-1} u_{j-1} + a_0 u_j + a_1 u_{j+1}. $$

In order to show that the norm of $A$ is less than 1, it is sufficient to show the following inequality:

$$ \sum_{j=0}^{J} v_j^2 \le \sum_{j=0}^{J} u_j^2. \qquad (11) $$

Furthermore, from the expression of $v_j$ and straightforward algebraic manipulations, we find the relation<fnref n="1" />:

$$ \begin{aligned} v_j^2 - u_j^2 = & - \frac{\nu - (\lambda a)^2}{2} \left( (u_j - u_{j-1})^2 + (u_{j+1} - u_j)^2 \right) + \frac{\nu^2 - (\lambda a)^2}{4} (u_{j+1} - 2u_j + u_{j-1})^2 \\ & - (\lambda a) u_j^2 + \frac{\nu(1 - \lambda a)}{2} (u_{j+1} - u_j)^2 + (\nu - \lambda a) u_j (u_{j+1} - u_j) \\ & + (\lambda a) u_{j-1}^2 - \frac{\nu(1 - \lambda a)}{2} (u_j - u_{j-1})^2 - (\nu - \lambda a) u_{j-1} (u_j - u_{j-1}), \end{aligned} $$

where the two last lines on the right-hand side are telescopic with respect to $j$. We sum this relation with respect to $j$ from 0 to $J$ and make use of the relations $u_{-1} = 0$, $u_{J+1} = u_J$ to simplify the boundary terms (that are obtained after summing the last two lines). We get:

$$ \begin{aligned} \sum_{j=0}^{J} v_j^2 - \sum_{j=0}^{J} u_j^2 = & - \frac{\nu - (\lambda a)^2}{2} \sum_{j=0}^{J} (u_j - u_{j-1})^2 + (u_{j+1} - u_j)^2 \\ & + \frac{\nu^2 - (\lambda a)^2}{4} \sum_{j=0}^{J} (u_{j+1} - 2u_j + u_{j-1})^2 - (\lambda a) u_J^2 - \frac{\nu(1 - \lambda a)}{2} u_0^2. \end{aligned} $$

<footnote n="1">¹More general decompositions of the same kind can be found in [3].</footnote>
From our assumption on $\lambda a$ and $\nu$, the very two last terms in the second line on the right-hand side are nonpositive, so we have:

$$ \sum_{j=0}^{J} v_j^2 - \sum_{j=0}^{J} u_j^2 \leq -\frac{\nu - (\lambda a)^2}{2} \sum_{j=0}^{J} (u_j - u_{j-1})^2 + (u_{j+1} - u_j)^2 \\ + \frac{\nu^2 - (\lambda a)^2}{4} \sum_{j=0}^{J} (u_{j+1} - 2u_j + u_{j-1})^2. $$

The inequality (11) follows in a straightforward way in the case $\nu \in [(\lambda a)^2, \lambda a]$ for in that case the right-hand side is the sum of two nonpositive quantities. We thus now examine the final case $\nu \in [\lambda a, 1]$, and we use the inequality<fnref n="2" />:

$$ (u_{j+1} - 2u_j + u_{j-1})^2 \leq 2(u_j - u_{j-1})^2 + 2(u_{j+1} - u_j)^2, $$

to get:

$$ \sum_{j=0}^{J} v_j^2 - \sum_{j=0}^{J} u_j^2 = -\frac{\nu - \nu^2}{2} \sum_{j=0}^{J} (u_j - u_{j-1})^2 + (u_{j+1} - u_j)^2 \leq 0, $$

which shows again the validity of (11). The proof of Lemma 1 is complete. □

One can wonder whether the power boundedness of $A$ is linked to our choice $k = 1$. If we had chosen $k = 2$ in (6), which corresponds to a second order extrapolation at the outflow boundary, the corresponding matrix $A$ would have read:

$$ A = \begin{pmatrix} a_0 & a_1 & 0 & \cdots & \cdots & \cdots & 0 \\ a_{-1} & a_0 & a_1 & \cdots & & & & \vdots \\ 0 & \cdots & \cdots & \cdots & \cdots & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \cdots & \cdots & \vdots \\ \vdots & & \ddots & a_{-1} & a_0 & a_1 & 0 & \vdots \\ \vdots & & & 0 & a_{-1} & a_0 & a_1 & \vdots \\ 0 & \cdots & \cdots & \cdots & 0 & a_{-1} - a_1 & a_0 + 2a_1 & \vdots \end{pmatrix} \in \mathcal{M}_{J+1}(\mathbb{R}), $$

and it is then a new problem to determine whether $A$ is power bounded (under the same restrictions on $\lambda a$ and $\nu$ or under more severe restrictions). It is actually shown in [4] that the above matrix $A$ is power bounded under the very same restrictions on $\lambda a$ and $\nu$ as in Lemma 1. The proof however relies on a suitable modification of the $\ell^2$ norm close to the outflow boundary, which is in the same spirit as the analysis in [8]. This gives, once again, a bound for the powers of $A$ that is even uniform with respect to the dimension $J$ but the above argument does not extend in a straightforward way. We do not know whether such energy arguments can be used for three point schemes and any a priori given value of $k$, but the analysis in [1] suggests that the corresponding matrix $A$ should be power bounded. This is left to further study.

We discuss below in Section 3 why the simple scenario of Lemma 1 does not extend to arbitrary stencils, but before that, we shall review some known facts on the numerical scheme (4), (5), (6) for arbitrary $r, p$ and $k$.

<footnote n="2">²This is a mere consequence of the inequality $(a-b)^2 \leq 2a^2 + 2b^2$.</footnote>
## 2.2 The general case

This section summarizes the main results of [3]. In addition to (3), we shall from now on assume that the coefficients $a_\ell$ satisfy the following (von Neumann) stability condition.

**Assumption 1.** There holds:

$$ \sup_{\theta \in \mathbb{R}} \left| \sum_{\ell=-r}^{p} a_{\ell} e^{i \ell \theta} \right| = 1. $$

Under Assumption 1, it is known that the convolution operator:

$$ \mathcal{T} : (u_j)_{j \in \mathbb{Z}} \in \ell^2(\mathbb{Z}; \mathbb{R}) \quad \mapsto \quad \left( \sum_{\ell=-r}^{p} a_\ell u_{j+\ell} \right)_{j \in \mathbb{Z}} \in \ell^2(\mathbb{Z}; \mathbb{R}) \quad (12) $$

has norm 1 and is therefore power bounded. Furthermore, the one-sided problem with Dirichlet boundary conditions on the left is also known to be a contraction. Namely, it is also known (see [2, 12]) that the solution to the numerical scheme:

$$ \left\{ \begin{array}{ll} u_j^{n+1} = \sum_{\ell=-r}^{p} a_\ell u_{j+\ell}^n, & j \in \mathbb{N}, \quad n \in \mathbb{N}, \\ u_\mu^n = 0, & \mu = -r, \dots, -1, \quad n \in \mathbb{N}, \end{array} \right. \quad (13) $$

satisfies:

$$ \forall n \in \mathbb{N}, \quad \sum_{j \in \mathbb{N}} (u_j^{n+1})^2 \le \sum_{j \in \mathbb{N}} (u_j^n)^2. $$

This contraction property holds because the time iteration in (13) can be written under the form $\pi \mathcal{T} \pi$ where $\pi$ denotes the orthogonal projection in $\ell^2(\mathbb{Z}; \mathbb{R})$ on the sequences that are supported on $\mathbb{N}$, and $\mathcal{T}$ is the pure convolution operator given in (12) whose norm equals 1. We should therefore always keep in mind that prescribing the Dirichlet boundary condition makes the $\ell^2$ norm decrease since it acts as an orthogonal projection and therefore preserves stability.

On the other hand, the outflow problem with extrapolation boundary condition:

$$ \left\{ \begin{array}{ll} u_j^{n+1} = \sum_{\ell=-r}^{p} a_\ell u_{j+\ell}^n, & j \le J, \quad n \in \mathbb{N}, \\ D_{-}^k u_{j+\mu}^n = 0, & \mu = 1, \dots, p, \quad n \in \mathbb{N}, \end{array} \right. \quad (14) $$

has also been analyzed in [3] where it is shown (under Assumption 1 and the consistency conditions (3)) that there exists a constant $C$ (that only depends on the given extrapolation order $k$) such that, independently of the initial condition, there holds:

$$ \forall n \in \mathbb{N}, \quad \sum_{j \le J} (u_j^n)^2 \le C \sum_{j \le J} (u_j^0)^2, $$

for any solution to (14). In other words, the time evolution operator in (14) is power bounded on $\ell^2((-\infty, J]; \mathbb{R})$, even though it is not necessarily a contraction for the $\ell^2$ norm.

In view of the above two results, it is only the coupling between the Dirichlet and the Neumann condition that may create a large time exponential instability. Namely, the combination of the above two stability results (the one for Dirichlet and the one for the extrapolation condition) implies the following estimate for the matrix $A$ in (7) (see the main result in [3]):

$$ \forall n \in \mathbb{N}, \quad \|A^n\| \le C_0 e^{n C_0 / J}, $$
where the constant $C_0$ is independent of $n$ and $J$, and $\| \cdot \|$ is the matrix norm associated with the $\ell^2$ norm on vectors. In particular, the spectral radius of $A$ satisfies a bound of the form (here and from now on $\rho$ stands for the spectral radius of a matrix):

$$ \rho(A) \le 1 + \frac{C}{J}, $$

where $C$ is a constant that does not depend on $J$. Hence, if $A$ admits an unstable eigenvalue of modulus larger than 1, it can not be “too large”. Our goal is to make explicit whether such unstable eigenvalues do indeed occur.

# 3 Examples of large time amplification

## 3.1 The basic principles

We make use of the theory of oscillating wave packets and shall use many times below the concept of *group velocity*. This notion is discussed in details in [9] in the context of finite difference schemes and we refer to that reference for a detailed presentation. A convenient reference for the same notion in the context of partial differential equations is [11]. The application of these notions to numerical boundary conditions is the purpose of [10] and we shall also use the analysis of [10] as a black box in our (formal) arguments below. Namely, from Assumption 1, we know that there are some values of $\theta \in \mathbb{R}$ such that the so-called *amplification factor*:

$$ \sum_{\ell=-r}^{p} a_{\ell} e^{i\ell\theta} $$

has modulus 1. For instance, the consistency conditions (3) show that this is the case at $\theta = 0$. From the von Neumann condition (that is, Assumption 1), such modes are the only ones that are not exponentially damped by the numerical scheme. Let $\underline{\theta}$ be such a real number and let us set:

$$ \underline{z} := \sum_{\ell=-r}^{p} a_{\ell} e^{i\ell\underline{\theta}} \in \mathbb{S}^{1}. $$

Then the plane wave:

$$ (n, j) \in \mathbb{N} \times \mathbb{Z} \mapsto \underline{z}^{n} e^{i j \underline{\theta}}, $$

is a solution to the numerical scheme (4), at least far from the boundaries. Writing $j\underline{\theta} = (j\Delta x)\underline{\theta}/\Delta x$, we can view the mesh width $\Delta x$ as a small wavelength parameter. The theory developed in [9] shows that slowly modulated highly oscillating wave packets of the form:

$$ \Phi(j \Delta x - \mathbf{v}_g n \Delta t) \underline{z}^n e^{i j \underline{\theta}} \qquad (15) $$

may be propagated by the numerical scheme (4) for a well-chosen group velocity $\mathbf{v}_g$. The (smooth) function $\Phi$ describes the *envelope* of the oscillations. In the trivial case $\underline{\theta} = 0$, $\underline{z} = 1$, one recovers the propagation of a smooth profile by a stable and consistent approximation of the transport equation; the group velocity $\mathbf{v}_g$ equals the transport velocity $a$ of (1) in that case. The above expression (15) does not exactly provide with a solution to (4) but it gives its *leading behavior* if one has in mind that the solution to (4) has an asymptotic expansion with respect to the small parameter $\Delta x$. This is exactly the same kind of arguments as in the so-called WKB analysis of linear geometric optics [7].

The boundary conditions on either side of the interval, meaning either (5) or (6), will now give rise to wave packet *reflection* as follows:
* We start at time $t = 0$ with some slowly modulated highly oscillating wave packet (15) that is supported in the middle of the computation interval $(0, L)$ (take for instance an envelope function $\Phi$ that is supported in $[L/3, 2L/3]$).

* Assume that the numerical wave packet (15) has a positive group velocity $v_g$, so that it will eventually reach the right boundary $x = L$ at some positive time.

* When hitting the boundary, the wave packet acts as a forcing term in the Neumann condition (6) of the form $\alpha_\mu \underline{z}^n$, $\mu = 1, \dots, p$, where $\alpha_\mu$ is computed by applying the Neumann condition to the wave packet itself. One should then determine whether the numerical scheme (4) supports wave packets of the form:

$$ \underline{z}^n \kappa^j, $$

with $|\kappa| \ge 1$ that will be reflected “backwards”. In the limit case $|\kappa| = 1$, the associated group velocity of the wave packet should be nonpositive so that the wave packet will propagate towards the left.

* In the case where in the previous reflection process, one of the group velocities is negative, the corresponding reflected wave packet will eventually hit the left boundary $x = 0$ and a similar reflection process will take place (the corresponding selection for the $\kappa$'s is now $|\kappa| \le 1$, and in the limit case $|\kappa| = 1$, the associated group velocity of the wave packet should be nonnegative).

Several points should be kept in mind: if a reflection gives rise to several wave packets, they will each carry part of the energy of the solution. In order to make the $\ell^2$ norm of the whole solution increase, it is desirable to have a sequence of reflections such that after one reflection on the right and one on the left, a wave packet with a given frequency $\theta$ and positive group velocity has its amplitude multiplied by a factor that has modulus larger than 1. Unless some cancelations appear, the wave packet (15) will necessarily be generated after reflection on the left of the interval since it has a positive group velocity. Observe now that it takes $O(J)$ time steps to make the whole transport back and forth from left to right, so an amplification pattern as described above should lead to a solution whose $\ell^2$ norm grows at least like $\exp(c n/J)$ for some positive constant $c$. This is precisely the situation which we try to highlight here for, in that case, we expect the spectral radius of $A$ to be larger than 1.

In what follows, we seek for such exponential amplifications by using as an initial condition a smooth profile, that is $\theta = 0$, $\underline{z} = 1$. To make the above mechanism work, the numerical scheme should admit an oscillating wave packet:

$$ 1^n \kappa^j $$

with $|\kappa| = 1$ and a negative group velocity ($\kappa$ is necessarily not equal to 1 for in that case the group velocity equals $a$ and is therefore positive). However, one should observe that when the Neumann condition (6) is applied to a smooth wave of the form $\Phi(j \Delta x - a n \Delta t)$, the corresponding error is $O(\Delta x^k)$ so the forcing term in (6) will have a very small reflection coefficient. Since the Dirichlet condition makes the $\ell^2$ norm decrease, this will not be sufficient to trigger an instability. We therefore need another auxiliary wave packet that is associated with a nonzero frequency and that has a positive group velocity.

Let us therefore summarize what we need to display some two boundary large time instability:

1. A finite difference scheme (4) that satisfies the conditions (3) and the von Neumann stability condition (Assumption 1).
<figure><img src="image_3.png" /><figcaption>Figure 1: Left: amplification factor (dark red curve) $\theta \mapsto \sum_{\ell=-r}^{p} a_{\ell} e^{i\ell\theta}$ with $r = p = 7$ and coefficients $a_{\ell}$ given in (16) within the unit circle $S^1$ (blue curve). Right: zoom of the amplification factor near the tangency point at $z = 1$.</figcaption></figure>

2. The finite difference scheme should also admit one plane wave of the form $1^n \kappa_1^j$ that is associated with a negative group velocity $v_{g,1}$ and another plane wave of the form $1^n \kappa_2^j$ that is associated with a positive group velocity $v_{g,2}$, both $\kappa_1$ and $\kappa_2$ being distinct from 1.

3. If the reflection coefficient from the first wave to the second on the right boundary is denoted $\alpha_{1\to2}$, and if the reflection coefficient from the second wave to the first on the right boundary is denoted $\beta_{2\to1}$, then we wish $|\alpha_{1\to2}\beta_{2\to1}| > 1$.

## 3.2 An example with first order extrapolation

In this first example, we have $r = p = 7$. The scheme reads as in (4) with the following choice for the coefficients:

$$a_{-7} = \frac{20133}{704759}, \qquad a_{-6} = -\frac{30433}{894654}, \qquad a_{-5} = \frac{20476}{371655}, \qquad a_{-4} = \frac{12703}{838076}, \qquad (16a)$$

$$a_{-3} = \frac{75599}{770583}, \qquad a_{-2} = \frac{31384}{945409}, \qquad a_{-1} = -\frac{4015}{85641}, \qquad a_{0} = \frac{261251}{274289}, \qquad (16b)$$

$$a_1 = \frac{53837}{928392}, \qquad a_2 = -\frac{27478}{587215}, \qquad a_3 = -\frac{54064}{714213}, \qquad a_4 = -\frac{27635}{674698}, \qquad (16c)$$

$$a_5 = -\frac{31244}{798847}, \qquad a_6 = \frac{23091}{711760}, \qquad a_7 = \frac{1864}{178339}. \qquad (16d)$$

The corresponding amplification factor is depicted in Figure 1. We implement this numerical scheme with the Dirichlet boundary condition (5) on the left of the interval and the Neumann condition (6) on the right ($k = 1$), that is:

$$\forall n \in \mathbb{N}, \quad u_{j+1}^n = \cdots = u_{j+7}^n = u_j^n.$$

The above numerical scheme is devised in such a way that it supports the following wave packets with
<figure><img src="image_4.png" /><figcaption>Figure 2: Evolution of $\ln \|\mathbf{u}^n\|_{\ell^2}$ as a function of $n\Delta t$ for the solution $\mathbf{u}^n = (u_j^n)_{j=0,\dots,J}$ of the numerical scheme starting from a smooth initial condition with the Dirichlet boundary condition (5) on the left of the interval and the Neumann condition (6) on the right ($k=1$) together with the coefficients given in (16). Here, we have set $J = 994$ with $\Delta t = \Delta x = 1/(J+1)$.</figcaption></figure>

given frequencies and group velocities:

$$ \begin{align*} (\underline{\kappa}_1, \mathbf{v}_{g,1}) &= (1, 1), \\ (\underline{\kappa}_2, \mathbf{v}_{g,2}) &= (e^{i0.288\pi}, -1), & (\underline{\kappa}_3, \mathbf{v}_{g,3}) &= (e^{-i0.288\pi}, -1), \\ (\underline{\kappa}_4, \mathbf{v}_{g,4}) &= (e^{i0.82\pi}, 1), & (\underline{\kappa}_5, \mathbf{v}_{g,5}) &= (e^{-i0.82\pi}, 1). \end{align*} $$

For all the above five values of $\kappa$, the associated amplification factor $z$ equals 1 so all wave packets are associated with the same time frequency, which allows for boundary reflection back and forth. If one starts from a smooth initial condition, the wave packet associated with the zero frequency $(\underline{\kappa}_1, \mathbf{v}_{g,1})$ will first travel to the right. Its reflection on the right boundary of the interval will give rise to the two wave packets associated with $(\underline{\kappa}_2, \mathbf{v}_{g,2})$ and $(\underline{\kappa}_3, \mathbf{v}_{g,3})$, both having small amplitudes since any smooth wave satisfies the Neumann condition with consistency error $O(\Delta x)$. When those reflected wave packets will hit the left boundary of the interval, the wave packet associated with $(\underline{\kappa}_1, \mathbf{v}_{g,1})$ will be generated again together with the wave packet associated with $(\underline{\kappa}_4, \mathbf{v}_{g,4})$ and $(\underline{\kappa}_5, \mathbf{v}_{g,5})$.

The amplification here will arise from the two pairs $(\underline{\kappa}_2, \underline{\kappa}_3)$ and $(\underline{\kappa}_4, \underline{\kappa}_5)$ of oscillating wave packets propagating back and forth from 0 to 1. We emphasize that the amplification is rather low, and we refer to Figure 2 for an illustration of the evolution of the logarithmic quantity $\ln \|\mathbf{u}^n\|_{\ell^2}$ for the solution of the numerical scheme $\mathbf{u}^n = (u_j^n)_{j=0,\dots,J}$ as a function of $n\Delta t$ when the numerical scheme is initialized with the following sequence:

$$ \forall j = 0, \dots, J, \quad u_j^0 = f(x_j), \quad \text{with } f(x) := \exp \left( -50 \left( x - \frac{1}{2} \right)^2 \right). $$

Here, we have used the following $\ell^2$ norm:

$$ \|\mathbf{u}^n\|_{\ell^2} := \left( \Delta x \sum_{j=0}^{J} |u_j^n|^2 \right)^{1/2} . $$
<figure><img src="image_5.png" /><figcaption>Figure 3: Left: space-time evolution of $|u_j^n|$, solution of the numerical scheme with the Dirichlet boundary condition (5) on the left of the interval and the Neumann condition (6) on the right ($k = 1$) together with the coefficients given in (16). Right: evolution of $\ln \|\mathbf{u}^n\|_{\ell^2}$ as a function of $n\Delta t$. Here, we have set $J = 994$ with $\Delta t = \Delta x = 1/(J+1)$.</figcaption></figure>

We observe, as expected, that it takes $O(1/\Delta x)$ amount of time for the instability to take off. The numerically computed slope of the linear regression of $\ln \|\mathbf{u}^n\|_{\ell^2}$ as a function of $n\Delta t$ is $1.63818 \times 10^{-2}$ which compares very well to:

$$ \frac{\rho(A) - 1}{\Delta x} \sim 1.63995 \times 10^{-2}, $$

where we denoted by $\rho(A)$ the spectral radius of the matrix $A$ of the numerical scheme. This phenomenon can be expected from the so-called power method for computing the largest eigenvalue of a matrix.

In order to better illustrate the reflection of the wave packets, Figure 3 shows the evolution of the solution of the numerical scheme initialized with the following (slowly modulated, highly oscillating) sequence:

$$ \forall j = 0, \dots, J, \quad u_j^0 = f(x_j), \quad \text{with } f(x) := \cos\left(\frac{0.82\pi}{\Delta x}\left(x - \frac{1}{2}\right)\right) \exp\left(-50\left(x - \frac{1}{2}\right)^2\right), $$

which corresponds to the linear superposition of the slowly modulated highly oscillating wave packets ($\kappa_4, \mathbf{v}_{g,4}$) and ($\kappa_5, \mathbf{v}_{g,5}$). As expected, we observe a reflection of the wave packets with an amplification over time.
<figure><img src="image_6.png" /><figcaption>Figure 4: Left: amplification factor (dark red curve) $\theta \mapsto \sum_{\ell=-r}^{p} a_{\ell} e^{i\ell\theta}$ with $r = p = 7$ and coefficients $a_{\ell}$ given in (17) within the unit circle $S^1$ (blue curve). Right: zoom of the amplification factor near the tangency point at $z = 1$.</figcaption></figure>

### 3.3 An example with second order extrapolation

In this second example, we still have $r = p = 7$. The scheme reads as in (4) with the following choice for the coefficients:

$$a_{-7} = \frac{17151}{869039}, \quad a_{-6} = -\frac{9591}{473236}, \quad a_{-5} = \frac{55269}{926798}, \quad a_{-4} = -\frac{1854}{92119}, \quad (17a)$$

$$a_{-3} = \frac{8983}{71254}, \quad a_{-2} = \frac{32579}{620922}, \quad a_{-1} = -\frac{47474}{813983}, \quad a_{0} = \frac{739673}{796040}, \quad (17b)$$

$$a_1 = \frac{2966}{34841}, \quad a_2 = -\frac{19830}{397889}, \quad a_3 = -\frac{71152}{650153}, \quad a_4 = -\frac{21338}{716071}, \quad (17c)$$

$$a_5 = -\frac{10189}{548431}, \quad a_6 = \frac{19029}{761263}, \quad a_7 = \frac{8820}{964529}. \quad (17d)$$

The corresponding amplification factor is depicted in Figure 4. We implement this numerical scheme with the Dirichlet boundary condition (5) on the left of the interval and the second order extrapolation condition (6) on the right ($k = 2$), that is:

$$\forall n \in \mathbb{N}, \quad u_{j+1}^n = 2u_j^n - u_{j-1}^n, \quad u_{j+2}^n = 2u_{j+1}^n - u_j^n \dots$$

which can be further simplified into:

$$\forall n \in \mathbb{N}, \quad \forall \mu = 1, \dots, p, \quad u_{j+\mu}^n = (\mu + 1) u_j^n - \mu u_{j-1}^n.$$

The above numerical scheme is devised in such a way that it supports the following wave packets with given frequencies and group velocities:

$$(\underline{\kappa}_1, \mathbf{v}_{g,1}) = (1, 1),$$

$$(\underline{\kappa}_2, \mathbf{v}_{g,2}) = (e^{i0.3\pi}, -1),$$

$$(\underline{\kappa}_3, \mathbf{v}_{g,3}) = (e^{-i0.3\pi}, -1),$$

$$(\underline{\kappa}_4, \mathbf{v}_{g,4}) = (e^{i0.8\pi}, 1),$$

$$(\underline{\kappa}_5, \mathbf{v}_{g,5}) = (e^{-i0.8\pi}, 1).$$

For all the above five values of $\kappa$, the associated amplification factor $z$ equals 1. The wave packet generation will thus be entirely similar as in our first example. Only the numerical values are different.
<figure><img src="image_7.png" /><figcaption>Figure 5: Left: space-time evolution of $|u_j^n|$, solution of the numerical scheme with the Dirichlet boundary condition (5) on the left of the interval and the Neumann condition (6) on the right ($k=2$) together with the coefficients given in (17). Right: evolution of $\ln \|\mathbf{u}^n\|_{\ell^2}$ as a function of $n\Delta t$. Here, we have set $J = 1000$ with $\Delta t = \Delta x = 1/(J+1)$.</figcaption></figure>

The amplification here will arise from the two pairs of oscillating wave packets propagating back and forth from 0 to 1. We emphasize that the amplification is stronger here than for the first order extrapolation. As could be expected, the reason is that the spectral radius of A will be larger.

Figure 5 shows the evolution of the solution of the numerical scheme initialized with the following (slowly modulated, highly oscillating) sequence:

$$ \forall j = 0, \dots, J, \quad u_j^0 = f(x_j), \quad \text{with } f(x) := \cos\left(\frac{0.8\pi}{\Delta x}\left(x - \frac{1}{2}\right)\right) \exp\left(-50\left(x - \frac{1}{2}\right)^2\right), $$

which corresponds to the linear superposition of the slowly modulated highly oscillating wave packets ($\kappa_4$, $v_{g,4}$) and ($\kappa_5$, $v_{g,5}$). We observe, as expected, a reflection of the wave packets with an amplification over time. The numerically computed slope of the linear regression of $\ln \|\mathbf{u}^n\|_{\ell^2}$ as a function of $n\Delta t$ is $3.15834 \times 10^{-1}$ which compares very well to the computed value:

$$ \frac{\rho(A) - 1}{\Delta x} \sim 3.15884 \times 10^{-1}, $$

where we denoted once again by $\rho(A)$ the spectral radius of the matrix A of the numerical scheme. The amplification is indeed much stronger in the present case.

## 3.4 Discussion

The above examples are meant to show that even though Dirichlet and Neumann numerical boundary conditions seem perfectly legitimate, large time stability issues may be an issue, even for one-step stable, explicit finite difference approximations of the transport equation. Still, the instability may be difficult to capture and it is quite unstable with respect to the various parameters, including the number of points $J$. If a large time instability occurs for some integer $J$, it may very well be that for $J+1$ the instability disappears. A connection with the spectral analysis of large Toeplitz matrices should be made.
Even though the reflection instability pattern is not hard to understand, there is no guarantee that any finite difference approximation with enough oscillating wave packets and appropriate group velocities will lead to an instability. To some extent, we have been lucky enough to obtain examples for $k = 1$ and $k = 2$.

Our main message is that the above examples show that the stability criterion exhibited in [1] may fail to be satisfied in several situations. This motivates setting numerical routines for trying to verify this criterion on high order schemes and more elaborate numerical boundary conditions.

## References

[1] Benoit, A. (2022). Stability of finite difference schemes approximation for hyperbolic boundary value problems in an interval. *Math. Comput.* **91**(335), 1171–1212.

[2] Coulombel, J.-F. and A. Gloria (2011). Semigroup stability of finite difference schemes for multidimensional hyperbolic initial-boundary value problems. *Math. Comput.* **80**(273), 165–203.

[3] Coulombel, J.-F. and F. Lagoutière (2020). The Neumann numerical boundary condition for transport equations. *Kinet. Relat. Models* **13**(1), 1–32.

[4] Coulombel, J.-F. and T. Lundquist (2020). On some stable boundary closures of finite difference schemes for the transport equation. Preprint, arXiv:2001.02583 [math.NA] (2020).

[5] Goldberg, M. (1977). On a boundary extrapolation theorem by Kreiss. *Math. Comput.* **31**, 469–477.

[6] Kreiss, H.-O. (1966). Difference approximations for hyperbolic differential equations. Numer. Solution partial diff. Equations, Proc. Sympos. Univ. Maryland 1965, 51-58 (1966).

[7] Lax, P. D. (1957). Asymptotic solutions of oscillatory initial value problems. *Duke Math. J.* **24**, 627–646.

[8] Strand, B. (1994). Summation by parts for finite difference approximations for $d/dx$. *J. Comput. Phys.* **110**(1), 47–67.

[9] Trefethen, L. N. (1982). Group velocity in finite difference schemes. *SIAM Rev.* **24**, 113–136.

[10] Trefethen, L. N. (1983). Group velocity interpretation of the stability theory of Gustafsson, Kreiss, and Sundstroem. *J. Comput. Phys.* **49**, 199–217.

[11] Whitham, G. B. (1999). *Linear and nonlinear waves*. (Paperback ed. ed.). Pure Appl. Math., Wiley-Intersci. Ser. Texts Monogr. Tracts. New York, NY: Wiley.

[12] Wu, L. (1995). The semigroup stability of the difference approximations for initial-boundary value problems. *Math. Comput.* **64**(209), 71–88.