# The Fourier spectral approach to the spatial discretization of quasilinear hyperbolic systems

Vincent Duchêne, Johanna Ulvedal Marstrander

► To cite this version:

Vincent Duchêne, Johanna Ulvedal Marstrander. The Fourier spectral approach to the spatial discretization of quasilinear hyperbolic systems. 2025. hal-05140303v2

HAL Id: hal-05140303

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

Preprint submitted on 5 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
# The Fourier spectral approach to the spatial discretization of quasilinear hyperbolic systems

Vincent Duchêne¹ and Johanna Ulvedal Marstrander²

¹IRMAR, Univ. Rennes, F-35000 Rennes, France.

²Department of Mathematical Sciences, NTNU, 7491 Trondheim, Norway

November 5, 2025

## Abstract

We discuss the rigorous justification of the spatial discretization by means of Fourier spectral methods of quasi-linear first-order hyperbolic systems. We provide uniform stability estimates that grant spectral convergence of the (spatially) semi-discretized solutions towards the corresponding continuous solution provided that the underlying system satisfies some suitable structural assumptions. We consider a setting with sharp low-pass filters and a setting with smooth low-pass filters and argue that —at least theoretically— smooth low-pass filters are operable on a larger class of systems. While our theoretical results are supported with numerical evidence, we also pinpoint some behavior of the numerical method that currently has no theoretical explanation.

**Keywords.** fourier spectral methods, spectral convergence, hyperbolic systems

**MSC codes.** 65M12, 65M70, 76M22

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>1.1</td><td>Motivation and related works</td><td>3</td></tr><tr><td>1.2</td><td>Definitions and notations</td><td>5</td></tr><tr><td>2</td><td>Symmetric quasilinear systems</td><td>6</td></tr><tr><td>2.1</td><td>Discretization with sharp low-pass filters</td><td>6</td></tr><tr><td>2.2</td><td>Discretization with smooth low-pass filters</td><td>10</td></tr><tr><td>3</td><td>Symmetrizable quasilinear systems</td><td>11</td></tr><tr><td>3.1</td><td>Discretization with smooth low-pass filters</td><td>12</td></tr><tr><td>3.2</td><td>Discretization with sharp low-pass filters</td><td>16</td></tr><tr><td>4</td><td>Numerical experiments for the Saint-Venant system</td><td>18</td></tr><tr><td>4.1</td><td>Analysis of the Saint-Venant system</td><td>19</td></tr><tr><td>4.2</td><td>Numerical experiments in dimension one</td><td>21</td></tr><tr><td>4.3</td><td>Numerical experiments in dimension two</td><td>24</td></tr><tr><td>A</td><td>Technical tools</td><td>26</td></tr></table>
# 1 Introduction

In this work we shall consider the spatial discretization by means of Fourier spectral methods of systems of the form

$$ \partial_t U + \sum_{j=1}^{d} A_j(U) \partial_{x_j} U = 0, \quad U|_{t=0} = U^0. \qquad (1.1) $$

where for all $j \in \{1, \dots, d\}$ and for all $U \in \mathbb{R}^n$, $A_j(U)$ are matrices satisfying Assumption A.1.

**Assumption A.1.** For all $j \in \{1, \dots, d\}$ and $U \in \mathbb{R}^n$, $A_j(U)$ are $n$-by-$n$ real-valued matrices. We assume that all entries of $A_j(\cdot)$ are polynomial.

**Remark 1.1.** Assuming that the system (1.1) has only polynomial nonlinearities may seem an over-restrictive assumption. This assumption is motivated by two considerations. Firstly, the Fourier spectral method may be efficiently implemented only within this framework; see Remark 1.2 below. Secondly, as far as we know, the literature lacks a theory analogous to para-differential calculus (see e.g. [19, Chapter 5]) on Sobolev spaces of periodic functions, which would provide the composition estimates in Proposition A.4 for general composition functions. All our results apply without assuming polynomial nonlinearities if Proposition A.4 holds without that assumption.

We shall assume further additional structural assumptions on the system (1.1), depending on the needs. In particular we shall always consider Friedrichs-symmetrizable systems, which guarantees the hyperbolicity of the system, and local-in-time well-posedness of the initial-value problem in $L^2$-based Sobolev spaces of sufficiently high regularity index, $U^0 \in H^s((2\pi\mathbb{T})^d)^n$ with $s > d/2 + 1$ (see Subsection 1.2). For simplicity, we consider in this work $2\pi$-periodic functions in all spatial directions. The results extend straightforwardly to more general periodic frameworks, and analogous results in the $n$-dimensional Euclidean space could be obtained with some simple adaptations.

Let $\mathcal{T}_N$ be the space of trigonometric polynomials of degree $N$:

$$ \mathcal{T}_N := \operatorname{span}\{\exp(i\boldsymbol{k} \cdot \boldsymbol{x}), \boldsymbol{k} \in \mathbb{Z}^d, |k_j| \le N, j = 1, \dots, d\}, \quad \text{and} \quad \mathcal{T}_N^n := \underbrace{\mathcal{T}_N \times \dots \times \mathcal{T}_N}_{n \text{ times}}. $$

Let $P_N: L^2((2\pi\mathbb{T})^d)^n \to \mathcal{T}_N^n$ be the $L^2$-projection operator onto $\mathcal{T}_N^n$: $P_N = \operatorname{Diag}(P_N(D))$, where $P_N$ is a Fourier multiplier with symbol $P_N(\cdot) = \mathbf{1}_{\lfloor -N, N \rfloor^d}(\cdot)$. We refer to $P_N$ as a **sharp low-pass filter**. For any $0 \le r \le s$ and $U \in H^s((2\pi\mathbb{T})^d)^n$ Sobolev space of order $s$, $P_N$ satisfies (see e.g. [5, (5.1.10),(5.8.4)])

$$ |U - P_N U|_{H^r} \le C(d, s, r) \langle N \rangle^{r-s} |U|_{H^s}. \qquad (1.2) $$

This kind of estimates are referred to in the literature as spectral convergence, and we will follow this terminology.

The standard Fourier spectral method (see e.g. [16, 13, 18] or [1, Section 3]) for the spatial discretization of the problem (1.1) amounts to seeking solutions $U_N: t \mapsto \mathcal{T}_N^n$ to the problem

$$ \partial_t U_N + P_N \left( \sum_{j=1}^{d} A_j(U_N) \partial_{x_j} U_N \right) = 0, \quad U_N|_{t=0} = P_N U^0. \qquad (1.3) $$

**Remark 1.2.** Let us recall that one of the great assets of Fourier spectral methods is that spatial differentiation and multiplication can be very efficiently performed (up to machine precision) by means of Fast Fourier Transform (FFT/IFFT) and multiplication at spatial collocations points, provided suitable dealiasing operations are performed. In practice, if the entries of $A_j(\cdot)$ are all polynomial with maximal degree $p$, then one computes (following Orszag's rule [20]) $N^{p+2}/2$ modes, and applying the projection $P_N$ after multiplications at collocation points performs the necessary dealiasing. If entries of $A_j(\cdot)$ are not polynomials, then one typically uses pseudo-spectral schemes that follow the aforementioned strategy but cannot be formulated as in (1.3).

We do not discuss in this work the full time-space discretization of (1.1), that is well-suited numerical time integrators for (1.3). The implementation we used for numerical experiments is described in more detail in Subsections 4.2 and 4.3.
In this work we shall discuss the long-time (i.e. uniform with respect to $N$) existence of solutions to (1.3) and spectral convergence towards solutions to (1.1) as $N \to \infty$.

Our results will depend on the structure of the system. As we shall see, the solution to (1.3) converges towards the corresponding solution to the problem (1.1) (assuming sufficient regularity) whenever the system is symmetric. If the system is only symmetrizable, the situation is more complicated. In order to deal with this situation we consider smooth low-pass filters, $S_N: L^2((2\pi\mathbb{T})^d)^n \to \mathcal{T}_N^n$ where $S_N = \text{Diag}(S_N(D))$ is a Fourier multiplier with symbol $S_N(\cdot) = S(\cdot/N)$ where $S$ is even and satisfies

$$ 
\begin{cases}
S(\boldsymbol{k}) = 1 & \text{if } \max_{j=1,...,d} |k_j| \le 1/2, \\
S(\boldsymbol{k}) = 0 & \text{if } \min_{j=1,...,d} |k_j| \ge 1, \\
S(\boldsymbol{k}) \in [0,1] & \text{otherwise,}
\end{cases}
$$

and $S^{1/2}$ is Lipschitz-continuous. When $d=1$, an example of such a function is $S_1(\cdot) := \max(0, \min(1, 2-2|\cdot|))^2$. When $d \ge 2$, one can set $S_d((k_1, \dots, k_d)) := S_1(k_1) \times \dots \times S_1(k_d)$. The advantage of such smooth low-pass filters is that —contrarily to the sharp low-pass filter— they satisfy commutator estimates with gains of regularity uniformly with respect to $N$; see Proposition A.6. These are crucial to parts of the analysis. Notice that because $|1-S(\boldsymbol{k})| \le 1 = |1-P_{1/2}(\boldsymbol{k})|$ when $\max_{j=1,...,d} |k_j| \ge 1/2$ and $1-S(\boldsymbol{k}) = 0 = 1-P_{1/2}(\boldsymbol{k})$ otherwise, we infer from (1.2) the corresponding spectral convergence estimate

$$ |U - S_N U|_{H^r} \le |U - P_{N/2} U|_{H^r} \le C(d, s, r) \langle N \rangle^{r-s} |U|_{H^s}. \quad (1.4) $$

The spatial discretization of (1.1) using smooth low-pass filters could amount to finding a solution $U_N: t \mapsto \mathcal{T}_N^n$ to

$$ \partial_t U_N + S_N \left( \sum_{j=1}^{d} A_j(U_N) \partial_{x_j} U_N \right) = 0, \quad U_N|_{t=0} = P_N U^0. \quad (1.5) $$

We also consider variants of this system such as

$$ \partial_t U_N + \sum_{j=1}^{d} (A_j^0 + S_N(A_j^1(U_N)[\circ])) \partial_{x_j} U_N = 0, \quad U_N|_{t=0} = P_N U^0. \quad (1.6) $$

where $A_j^0 = A_j(0)$ and $A_j^1(U) = A_j(U) - A_j^0$. Indeed, applying smooth low-pass filters to linear terms is unnecessary, especially when one uses exponential time integrators (see e.g. Program 27 in [25]). Notice the distinction between (1.5) and (1.6) is only necessary when using smooth low-pass filters since when using the sharp low-pass filter $P_N$ one has $P_N U_N = U_N$ and $P_N$ commutes with $A_j^0$.

**Outline** Let us now describe the structure and main results of this work. Symmetric systems are discussed in Section 2. We consider sharp and smooth low-pass filters in Subsection 2.1 and 2.2 respectively. We obtain convergence of the semi-discretized solutions in both cases, stated in Propositions 2.4 and 2.6. Symmetrizable systems are discussed in Section 3. Subsection 3.1 concerns smooth low-pass filters and we obtain analogous convergence results, stated in Proposition 3.5. The case with sharp low-pass filters is treated in Subsection 3.2 and in order to secure spectral convergence, more stringent structural assumptions on the system are required. This yields Proposition 3.10. Numerical experiments illustrating and investigating the sharpness of our theoretical results are provided in Section 4.

## 1.1 Motivation and related works

Our work was motivated by the study of Boussinesq and Whitham-Boussinesq systems that are nonlinear dispersive models for the propagation of surface gravity waves [17, 7]. The Fourier spectral method is especially indicated for the spatial discretization of these systems since the nonlinear contributions are quadratic (recall Remark 1.2) and the
dispersive contributions take the form of Fourier multipliers. The spectral convergence of discretized versions of Boussinesq models towards the corresponding continuous solutions was proved in [26] and [6]. These results have a shortcoming in that they lack uniformity in the non-dispersive limit. In order to clarify this point, let us consider the Benjamin–Bona–Mahony (BBM) and (a variant of) the Whitham equations, which read respectively

$$ (\mathrm{Id} - \mu \partial_x^2) \partial_t u + \partial_x u + u \partial_x u = 0, \quad \text{and} \quad \partial_t u + \frac{\tanh(\sqrt{\mu}|D|)}{\sqrt{\mu}|D|} (\partial_x u + u \partial_x u) = 0. $$

The aforementioned Boussinesq (resp. Whitham–Boussinesq) systems can be loosely considered as systems extending the BBM (resp. Whitham) scalar equation, in the same way the Saint-Venant system discussed in Section 4 extends the inviscid Burgers equation. All these equations provide valid approximations of water waves provided (among other assumptions) the shallowness parameter $\mu > 0$ satisfies $\mu \ll 1$; see [17, 7].

Results proved in [26] and [6] and adapted to the simplest case of scalar equations take the form

$$ |u_N(t, \cdot)|_{H^s} \le C_\mu(t) |u_N|_{t=0}|_{H^s} \quad \text{and} \quad |(u - u_N)(t, \cdot)|_{H^r} \le C_\mu(t) N^{r-s} \quad \text{for any } s > 1/2 \text{ and } 0 \le r \le s, $$

where $u$ is the solution to the scalar equation, and $u_N$ the corresponding solution to the semi-discretized equation

$$ (\mathrm{Id} - \mu \partial_x^2) \partial_t u_N + \partial_x u_N + P_N(u_N \partial_x u_N) = 0, \quad \text{and} \quad \partial_t u_N + \frac{\tanh(\sqrt{\mu}|D|)}{\sqrt{\mu}|D|} (\partial_x u_N + P_N(u_N \partial_x u_N)) = 0. $$

The aforementioned shortcoming is that $C_\mu(t)$ depends nonuniformly on $\mu$ in the non-dispersive limit $\mu \ll 1$, typically $C_\mu(t) \lesssim \exp(\mu^{-1}t)$. This is inconsistent with the standard well-posedness theory for initial data in Sobolev spaces $H^s(2\pi\mathbb{T})$, $s > 3/2$, which holds uniformly with respect to $\mu \in (0, 1]$ (see e.g. [15, Proposition 6] and [22, 10, 21] for the corresponding results on the Boussinesq and Whitham–Boussinesq systems). The reason for this discrepancy is that for the spectral convergence results, stability estimates on the BBM and Whitham equations are obtained by viewing them as *semilinear systems* which can be formulated using the Duhamel formula

$$ u(t, \cdot) = \int_0^t e^{-(t-\tau)L_\mu(D)} L_\mu(D) \left( \frac{u(\tau, \cdot)^2}{2} \right) d\tau $$

where $L_\mu(D) = \frac{\partial_x}{\mathrm{Id} - \mu \partial_x^2}$ for the BBM equation, and $L_\mu(D) = \frac{\tanh(\sqrt{\mu}|D|)}{\sqrt{\mu}|D|} \partial_x$ for the Whitham equation. Notice that in both cases $L_\mu(D) \in \mathcal{B}(L^2(2\pi\mathbb{T}); L^2(2\pi\mathbb{T}))$, but that $\sup\{|L_\mu v|_{L^2} : |v|_{L^2} = 1\}$ is not uniformly bounded with respect to $\mu \in (0, 1]$, which is the source of the issue when $\mu \ll 1$. On the other hand, one can view the BBM and Whitham equation as perturbations of the inviscid Burgers equation, with skew-symmetric dispersive terms that are inconspicuous for the energy method. This leads to stability estimates in $H^s(2\pi\mathbb{T})$, $s > 3/2$, that are uniform with respect to $\mu \ll 1$. Of course obtaining such results on dispersive systems requires a good understanding of the underlying (non-dispersive) quasilinear systems, which is the focus of the current work.

The rigorous analysis of semi-discretized or fully (space and time) discretized Fourier spectral schemes for semi-linear equations is a very rich and active topic, which is impossible to summarize within a few lines; let us simply mention [4, 3] which are particularly relevant as they specifically consider the BBM equation and pay attention to the non-dispersive limit (although together with vanishing nonlinearity, that is the long wave limit). In contrast, to the best of our knowledge, there are only a handful of works dedicated to the rigorous analysis of the Fourier method for *quasilinear systems*, culminating with the work of Bardos and Tadmor [1] (following [23, 24, 11]); see also the review paper [12]). Here the authors consider the inviscid Burgers equation, as well as the one-dimensional isentropic Euler equation and the incompressible Euler equations. The first two systems belong to the class of equations we study. Specifically, the inviscid Burgers equation has the symmetric structure we employ in Subsection 2.1, while the isentropic Euler equation in Lagrangian coordinates enjoys the Hamiltonian structure discussed in Subsection 3.2. The incompressible Euler equations

$$ \partial_t u + \mathbb{P}((u \cdot \nabla)u) = 0, \quad \mathbb{P}(u) = u, \quad \mathbb{P}(u) := u - \nabla \Delta^{-1} \nabla \cdot u $$
does not belong to the class of equations studied in this work due to the presence of the Leray projection operator $\mathbb{P}$ but it would not be difficult to extend our analysis (specifically concerning the symmetric situation) to this system. Let us mention that while our analysis is very similar, the estimates obtained in [1] are not as sharp as the ones obtained in our work due to different choices when performing stability estimates, as the former take the form (see Theorem 3.1 therein)

$$|(U - U_N)(t, \cdot)|_{L^2} \le C(t)N^{-s}|U|_{t=0}|_{H^s} + N^{\frac{3}{4}-\frac{s}{2}} \max_{\tau \le t} |U(\tau, \cdot)|_{H^s},$$

while our results do not feature the second contribution. Quite interestingly, the authors in [1] also prove the *emergence of spurious oscillations* of the semi-discretized solution after the critical time of shock formation for the continuous solution of the inviscid Burgers equation. It would be interesting to study this problem when smooth low-pass filters are used and compare with the spectral viscosity method described therein.

## 1.2 Definitions and notations

In this section, we introduce a few notations used throughout the work.

Let $L^2((2\pi\mathbb{T})^d)$ be the Lebesgue space of real-valued, square-integrable functions on the $2\pi$-periodic torus and

$$L^2((2\pi\mathbb{T})^d)^n = \underbrace{L^2((2\pi\mathbb{T})^d) \times \dots \times L^2((2\pi\mathbb{T})^d)}_{n \text{ times}}.$$

We endow $L^2((2\pi\mathbb{T})^d)$ with the standard Lebesgue norm denoted $|\cdot|_{L^2}$, and the corresponding inner-product is denoted $(\cdot, \cdot)|_{L^2}$. Similarly, we denote by $L^\infty((2\pi\mathbb{T})^d)$ the Lebesgue space of bounded functions and $W^{1,\infty}((2\pi\mathbb{T})^d)$ the space of Lipschitz continuous functions, endowed with their natural norms.

We use the notation $\langle \cdot \rangle = (1 + |\cdot|^2)^{1/2}$ and $\Lambda^s = (1 - \Delta)^{s/2}$, i.e. the Fourier multiplier with symbol $\langle \cdot \rangle^s$ (see e.g. [19] for Fourier multipliers). For real $s \ge 0$, we denote the $L^2$-based periodic Sobolev spaces by $H^s((2\pi\mathbb{T})^d)^n$:

$$H^s((2\pi\mathbb{T})^d)^n = \{U \in L^2((2\pi\mathbb{T})^d)^n, |U|_{H^s} < \infty\}, \text{ where} \\ |U|_{H^s}^2 = |\Lambda^s U|_{L^2}^2 = \sum_{k \in \mathbb{Z}^d} \langle k \rangle^{2s} |\hat{U}_k|^2$$

where $\hat{U}_k$ is the $k$-th fourier coefficient of $U$. For $X$ a Banach space and $I \subset \mathbb{R}$ an interval, the space of $X$-valued continuous functions on $I$ is denoted $C(I; X)$. Given $n \in \mathbb{N}$, the space of continuously $n$-th differentiable functions is denoted $C^n(I; X)$.

As mentioned previously, we denote $\mathcal{T}_N$ the space of trigonometric polynomials of degree $N$:

$$\mathcal{T}_N := \operatorname{span}\{\exp(i\mathbf{k} \cdot \mathbf{x}), \mathbf{k} = (k_1, \dots, k_d) \in \mathbb{Z}^d, |k_j| \le N, j = 1, \dots, d\}, \text{ and } \mathcal{T}_N^n := \underbrace{\mathcal{T}_N \times \dots \times \mathcal{T}_N}_{n \text{ times}}$$

and $P_N := \operatorname{Diag}(P_N(D))$ the Fourier multiplier with symbol $P_N := \mathbf{1}_{[-N,N]^d}$ is the $L^2$-projection operator onto $\mathcal{T}_N^n$. Here, $[-N, N] := \{-N, -N+1, \dots, N-1, N\}$.

We set $S_N = \operatorname{Diag}(S_N(D))$ a Fourier multiplier with symbol $S_N(\cdot) = S(\cdot/N)$ where $S$ is even and satisfies

$$\left\{ \begin{array}{ll} S(\mathbf{k}) = 1 & \text{if } \max_{j=1,...,d} |k_j| \le 1/2, \\ S(\mathbf{k}) = 0 & \text{if } \min_{j=1,...,d} |k_j| \ge 1, \\ S(\mathbf{k}) \in [0,1] & \text{otherwise,} \end{array} \right.$$

and $S^{1/2}$ is Lipschitz-continuous. Apart from these properties, the specific profile of the symbol $S$ is inconsequential.

We denote by $C(\lambda_1, \lambda_2, \dots)$ a positive “constant” depending on its parameters. Whenever such a parameter represents the norm of a function, $C$ depends non-decreasingly on said norm. Whenever the parameter is a subset of the Euclidean space, $C$ depends non-decreasingly on this parameter when set inclusion is used as (partial) ordering. Dependency on regularity indices $s \in \mathbb{R}$ or the dimension $d$ are omitted when it is unessential or clear from the context. Note that such constants $C$ will always be independent of the degree of the approximation $N$.
## 2 Symmetric quasilinear systems

In this section, we consider systems (1.1), where for all $j \in \{1, \dots, d\}$, $A_j(U)$ satisfies the Assumption A.1 and is additionally self-adjoint. We have the following standard result (see e.g. [2]).

**Proposition 2.1 (Well-posedness).** Let $s > 1 + d/2$, and $M > 0$. Suppose that for all $j \in \{1, \dots, d\}$, $A_j(\cdot)$ satisfies the Assumption A.1 and is self-adjoint. There exists $C > 0$ and $T > 0$ (depending only on $s$ and $M$) such that for every $U^0 \in H^s((2\pi\mathbb{T})^d)^n$ such that $|U^0|_{H^s} \le M$, there exists a unique $U \in C(I; H^s((2\pi\mathbb{T})^d)^n)$ maximal-in-time classical solution to (1.1), and moreover the open time interval $I \supset [0, T/|U^0|_{H^s}]$ and for all $0 \le t \le T/|U^0|_{H^s}$,

$$|U|_{H^s} \le |U^0|_{H^s} \exp(C|U^0|_{H^s}t) \le 2|U^0|_{H^s}.$$

### 2.1 Discretization with sharp low-pass filters

Recall that the spatial discretization of the problem (1.1) amounts to finding a solution $U_N: t \mapsto T_N^n$ to the problem

$$\partial_t U_N + P_N \left( \sum_{j=1}^{d} A_j(U_N) \partial_{x_j} U_N \right) = 0, \quad U_N|_{t=0} = P_N U^0, \qquad (2.1)$$

where $P_N = \text{Diag}(P_N(D))$, with $P_N(\cdot) = \mathbf{1}_{[-N,N]}(\cdot)$.

We want to show that the semi-discretized solutions $U_N$ to (2.1) converge as $N \to \infty$ towards $U$ the corresponding solution to the underlying system (1.1). To do so, we will first show that the semi-discretized solutions exist and are bounded on a time interval independent of $N$. This is Proposition 2.2. Then we use this bound to compare the semi-discretized and continuous solution on the interval of existence in Proposition 2.4. Finally, we refine this result by showing that if $N$ is large enough the existence of the semi-discretized solution and the estimate on the difference hold on any compact subset of the interval of existence of the solution to (1.1).

**Proposition 2.2 (Uniform estimates).** Let $s > d/2 + 1$ and $M > 0$. Suppose that for all $j \in \{1, \dots, d\}$, $A_j(\cdot)$ satisfies the Assumption A.1 and is self-adjoint. There exists $C > 0$ and $T > 0$ depending only on $s$ and $M$ such that for every $N \in \mathbb{N}$ and for every $U^0 \in H^s((2\pi\mathbb{T})^d)^n$ such that $|U^0|_{H^s} \le M$, there exists a unique $U_N \in C(I_N; H^s((2\pi\mathbb{T})^d)^n)$ maximal-in-time classical solution to (2.1) and $U_N|_{t=0} = P_N U^0$. The open time interval $I_N \supset [0, T/|U^0|_{H^s}]$ and for all $0 \le t \le T/|U^0|_{H^s}$,

$$|U_N|_{H^s} \le |U^0|_{H^s} \exp(C|U^0|_{H^s}t) \le 2|U^0|_{H^s}.$$

Moreover, for any $s' \ge s$, one has $U_N \in C^1(I_N; H^s((2\pi\mathbb{T})^d)^n)$ and for any $0 < T^* \in I_N$ and $M^* > 0$ such that $\sup_{t \in [0, T^*]} |U_N(t, \cdot)|_{H^s} \le M^*$ there exists $C^* > 0$ depending only on $s$, $s'$ and $M^*$ such that for all $0 \le t \le T^*$,

$$|U_N|_{H^s'} \le |U^0|_{H^s'} \exp(C^*M^*t).$$

The key ingredient to show Proposition 2.2 is the following apriori estimate.

**Lemma 2.3.** Let $s' \ge s > d/2 + 1$. Suppose that for all $j \in \{1, \dots, d\}$, $A_j(\cdot)$ satisfies the Assumption A.1 and is self-adjoint. Let $U_N \in C(I_N; H^s((2\pi\mathbb{T})^d)^n)$ be solution to (2.1) on $I_N \subset \mathbb{R}$ an open time interval. Then $U_N \in C^1(I_N; H^s((2\pi\mathbb{T})^d)^n)$ and for any $t \in I_N$, one has

$$\frac{d}{dt} |U_N|_{H^s'} \le C (|U_N|_{H^s}) |U_N|_{H^s} |U_N|_{H^s'}. \qquad (2.2)$$

*Proof.* Notice $U_N = P_N U_N$ and hence we have smoothness in space, $U_N \in C(I_N; H^\sigma((2\pi\mathbb{T})^d)^n)$ for any $\sigma \in \mathbb{R}$. We infer smoothness in time, $U_N \in C^1(I_N; H^\sigma((2\pi\mathbb{T})^d)^n)$, using eq. (2.1), product and composition estimates
in $H^s((2\pi\mathbb{T})^d)^n$—Propositions A.3 and A.4—and that $P_N : H^s((2\pi\mathbb{T})^d)^n \to H^\sigma((2\pi\mathbb{T})^d)^n$ is bounded. Denote $\dot{U}_N := \Lambda^{s'} U_N$. Using that $U_N$ satisfies the system (2.1), that $P_N$ is symmetric for the $L^2$ inner-product and $P_N \dot{U}_N = \dot{U}_N$, and finally that $A_j$ are self-adjoint for all $j \in \{1, \dots, d\}$ and integration by parts, we have

$$ 
\begin{align*}
\frac{1}{2} \frac{d}{dt} |\dot{U}_N|^2_{L^2} &= (\Lambda^{s'} \partial_t U_N, \dot{U}_N)_{L^2} = -(P_N(\sum_{j=1}^d [\Lambda^{s'}, A_j(U_N)]\partial_{x_j} U_N), \dot{U}_N)_{L^2} - (P_N(\sum_{j=1}^d A_j(U_N)\partial_{x_j} \dot{U}_N), \dot{U}_N)_{L^2} \\
&= -(\sum_{j=1}^d [\Lambda^{s'}, A_j(U_N)]\partial_{x_j} U_N, \dot{U}_N)_{L^2} - (\sum_{j=1}^d A_j(U_N)\partial_{x_j} \dot{U}_N, \dot{U}_N)_{L^2} \\
&= -(\sum_{j=1}^d [\Lambda^{s'}, A_j(U_N)]\partial_{x_j} U_N, \dot{U}_N)_{L^2} + \frac{1}{2}(\sum_{j=1}^d [\partial_{x_j}, A_j(U_N)]\dot{U}_N, \dot{U}_N)_{L^2}.
\end{align*}
$$

Using the commutator and composition estimates of Propositions A.5,A.4 with $s_0 = s - 1$ for the first contribution, the continuous Sobolev embedding $H^s((2\pi\mathbb{T})^d)^n \subset W^{1,\infty}((2\pi\mathbb{T})^d)^n$ (Proposition A.1) for the second as well the Cauchy–Schwarz inequality gives the bound

$$ \frac{1}{2} \frac{d}{dt} |\dot{U}_N|^2_{L^2} \le C(|U_N|_{H^s}) |U_N|_{H^s} |\dot{U}_N|^2_{L^2}, $$

which yields the desired inequality. □

Now we proceed to prove Proposition 2.2.

*Proof of Proposition 2.2.* For each $N > 0$, existence and uniqueness of a local-in-time solution $U_N$ to (2.1) follows from its formulation as a system of ODEs in the Banach space $H^s((2\pi\mathbb{T})^d)^n$ (using eq. (2.1), product and composition estimates in $H^{s-1}((2\pi\mathbb{T})^d)^n$—Propositions A.3 and A.4—and that $P_N : H^{s-1}((2\pi\mathbb{T})^d)^n \to H^s((2\pi\mathbb{T})^d)^n$ is bounded) and the Picard–Lindelöf theorem. We denote $I_N$ the maximal interval of existence, and notice as in the proof of Lemma 2.3 that $U_N \in C^1(I_N; H^{s'}((2\pi\mathbb{T})^d)^n)$ for any $s' \ge s$.

The second part of the proposition is an immediate consequence of Lemma 2.3 and Grönwall's inequality, using that $|U_N|_{t=0}|_{H^{s'}} = |P_N U^0|_{H^{s'}} \le |U^0|_{H^{s'}}$. To show the first part of the proposition, we use a standard continuity argument. Let

$$ \varphi_N : t \mapsto \sup_{t' \in [0, t]} |U_N(t', \cdot)|_{H^s} \quad \text{and} \quad J_N := \{t \in I_N \cap \mathbb{R}^+ : \varphi_N(t) \le 2|U^0|_{H^s}\}. $$

Since $U_N \in C(I_N; H^s((2\pi\mathbb{T})^d)^n)$, $\varphi_N \in C(I_N \cap \mathbb{R}^+; \mathbb{R})$ is non-decreasing and $J_N = \varphi_N^{-1}([0, 2|U^0|_{H^s}])$ is a closed interval. Let us prove that one can set $T > 0$ independently of $N$ such that $J_N \cap [0, T/|U^0|_{H^s}]$ is a non-empty open subset of $[0, T/|U^0|_{H^s}]$. Notice $0 \in J_N$. Let $t \in J_N \cap [0, T/|U^0|_{H^s}]$. By the second part of the proposition, we can set $C = 2C^*$ depending only on $s$ and $|U^0|_{H^s}$ such that

$$ |U_N(t, \cdot)|_{H^s} \le |U^0|_{H^s} \exp(C|U^0|_{H^s} t) \le |U^0|_{H^s} \exp(CT). \quad (2.3) $$

Choosing $T = \ln(3/2)/C$, we find that $\varphi_N(t) \le \frac{3}{2}|U^0|_{H^s}$, and hence —by the continuity and monotonicity of $\varphi_N$— there exists $\delta > 0$ such that $[0, t + \delta] \in J_N$. This concludes the proof that $J_N \cap [0, T/|U^0|_{H^s}]$ is a non-empty, connected, closed and open subset of $[0, T/|U^0|_{H^s}]$. Hence $I_N \supset J_N \supset [0, T/|U^0|_{H^s}]$ and estimate (2.3) concludes the proof. □

Having established a bound on the semi-discretized solution $U_N$, we proceed to estimate the difference between the solution $U_N$ to (2.1) and $U$, solution to (1.1).
**Proposition 2.4 (Convergence).** Let $s > d/2 + 1$, $M > 0$. Suppose that $U^0 \in H^s((2\pi\mathbb{T})^d)^n$ with $|U^0|_{H^s} \le M$. Suppose that for all $j \in \{1, \dots, d\}$, $A_j(\cdot)$ satisfies the Assumption A.1 and is self-adjoint. Denote $U \in C(I; H^s((2\pi\mathbb{T})^d)^n)$ the maximal-in-time classical solution to (1.1), and $U_N \in C(I_N; H^s((2\pi\mathbb{T})^d)^n)$ the maximal-in-time classical solution to (2.1). Let $T > 0$ be the minimum value of Propositions 2.1 and 2.2. For all $0 \le t \le T/|U^0|_{H^s}$, there is a $C > 0$, depending only on $s$ and $M$ such that for any $0 \le r \le s$,

$$|(U - U_N)(t, \cdot)|_{H^r} \le C |U^0|_{H^s} N^{r-s}.$$

Moreover, for every compact subset $I^* \subset I$, there is an $N_0 \in \mathbb{N}$ and $C^* > 0$, depending only on $s$, $|I^*|$ and $M^* := \sup_{t \in I^*} |U(t, \cdot)|_{H^s}$ such that for all $N \ge N_0$, $I_N \supset I^*$ and for any $0 \le r \le s$,

$$\sup_{t \in I^*} |(U - U_N)(t, \cdot)|_{H^r} \le C^* M^* N^{r-s}.$$

*Proof.* We shall first prove the result for $t \in [0, T/|U^0|_{H^s}]$. Let us assume first that $U^0 \in H^{s+1}((2\pi\mathbb{T})^d)^n$ so that $U \in C(I; H^{s+1}((2\pi\mathbb{T})^d)^n) \cap C^1(I; H^s((2\pi\mathbb{T})^d)^n)$, the general case being deduced afterwards. Recall that $T > 0$ is the minimum value of Propositions 2.1 and 2.2, so that we have $|U(t, \cdot)|_{H^s} \le 2|U^0|_{H^s}$ and $|U_N(t, \cdot)|_{H^s} \le 2|U^0|_{H^s}$ for all $N \in \mathbb{N}$. Denote $D_N := U - U_N$ and notice

$$\partial_t D_N + \sum_{j=1}^{d} A_j(U_N) \partial_{x_j} D_N + \sum_{j=1}^{d} (A_j(U) - A_j(U_N)) \partial_{x_j} U = -(\mathrm{Id} - P_N) \left( \sum_{j=1}^{d} A_j(U_N) \partial_{x_j} U_N \right).$$

Now we apply the smooth low-pass filter $S_N$ and use that $S_N(\mathrm{Id} - P_N) = 0$:

$$\partial_t S_N D_N + S_N^{1/2} \left( \sum_{j=1}^d A_j(U_N) \partial_{x_j} S_N^{1/2} D_N \right) + S_N^{1/2} \left( \sum_{j=1}^d [S_N^{1/2}, A_j(U_N)] \partial_{x_j} D_N \right) \\ + S_N \left( \sum_{j=1}^d (A_j(U) - A_j(U_N)) \partial_{x_j} U \right) = 0.$$

Testing against $D_N$, using that $S_N$ is symmetric for the $L^2((2\pi\mathbb{T})^d)^n$ inner-product and commutes with $\partial_t$ and after some rearranging we get

$$\frac{d}{dt} |S_N^{1/2} D_N|_{L^2}^2 = - \left( \sum_{j=1}^d A_j(U_N) \partial_{x_j} S_N^{1/2} D_N, S_N^{1/2} D_N \right)_{L^2} - \left( \sum_{j=1}^d [S_N^{1/2}, A_j(U_N)] \partial_{x_j} D_N, S_N^{1/2} D_N \right)_{L^2} \\ - \left( S_N^{1/2} \left( \sum_{j=1}^d (A_j(U) - A_j(U_N)) \partial_{x_j} U \right), S_N^{1/2} D_N \right)_{L^2}.$$

Using integration by parts and the properties of $A_j(\cdot)$—Assumption A.1 together with the composition estimate of Proposition A.4 and self-adjointness—and the continuous Sobolev embedding $H^s((2\pi\mathbb{T})^d)^n \subset W^{1,\infty}((2\pi\mathbb{T})^d)^n$ (Proposition A.1), the first term on the right-hand side is bounded by $C(|U_N|_{H^s}) |U_N|_{H^s} |S_N^{1/2} D_N|_{L^2}^2$. Using the commutator estimate of Proposition A.6 and the composition estimate of Proposition A.4, the second term on the right is bounded by $C(|U_N|_{H^s}) |U_N|_{H^s} |D_N|_{L^2} |S_N^{1/2} D_N|_{L^2}$. By the composition estimate of Proposition A.4, the third term is bounded by $C(|U|_{H^s}, |U_N|_{H^s}) |U|_{H^s} |D_N|_{L^2} |S_N^{1/2} D_N|_{L^2}$, where we again use Assumption A.1 as well as the boundedness of $S_N^{1/2}$. In all three estimates, we also used the Cauchy–Schwarz inequality. Altogether, we get

$$\frac{1}{2} \frac{d}{dt} (|S_N^{1/2} D_N|_{L^2}^2) \le C(|U|_{H^s}, |U_N|_{H^s}) (|U_N|_{H^s} + |U|_{H^s}) (|S_N^{1/2} D_N|_{L^2}^2 + |D_N|_{L^2}^2) |S_N^{1/2} D_N|_{L^2}.$$
Now we remark that since the symbol $S_N$ satisfies $S_N(\cdot) \in [0, 1]$ and $S_N(k) = 1$ if $\max_{j=1,...,d} |k_j| \le N/2$ we have
$|(1 - S_N^{1/2}(\cdot))(\cdot)^{-s}|_{L^\infty} \le \langle N/2 \rangle^{-s}$ and hence

$$
|\mathbf{D}_N|_{L^2} \le |\mathbf{S}_N^{1/2} \mathbf{D}_N|_{L^2} + |(\mathrm{Id} - \mathbf{S}_N^{1/2}) \mathbf{D}_N|_{L^2} \le |\mathbf{S}_N^{1/2} \mathbf{D}_N|_{L^2} + \langle N/2 \rangle^{-s} |\mathbf{D}_N|_{H^s}.
$$

Hence since $|\mathbf{D}_N|_{H^s} \le |\mathbf{U}|_{H^s} + |\mathbf{U}_N|_{H^s} \le 4|\mathbf{U}^0|_{H^s}$ by the triangle inequality we have

$$
\frac{1}{2} \frac{d}{dt} (|S_N^{1/2} D_N|_{L^2}^2) \le C(|U^0|_{H^s}) |U^0|_{H^s} (|S_N^{1/2} D_N|_{L^2} + |U^0|_{H^s} N^{-s}) |S_N^{1/2} D_N|_{L^2}
$$

and we infer by Grönwall's Lemma that

$$
|S_N^{1/2} \mathbf{D}_N(t, \cdot)|_{L^2} \le C(|\mathbf{U}^0|_{H^s}) |\mathbf{U}^0|_{H^s} (|S_N^{1/2} \mathbf{D}_N|_{t=0}|_{L^2} + tN^{-s} |\mathbf{U}^0|_{H^s}) \exp \left( C(|\mathbf{U}^0|_{H^s}) |\mathbf{U}^0|_{H^s} t \right).
$$

The desired estimate for $r = 0$ follows by using that $S_N^{1/2} \mathbf{D}_N|_{t=0} = S_N^{1/2}(\mathrm{Id} - P_N)\mathbf{U}^0 = 0$, $t \in [0, T/|\mathbf{U}^0|_{H^s}]$ and again

$$
|\mathbf{D}_N|_{L^2} \le |S_N^{1/2} \mathbf{D}_N|_{L^2} + \langle N/2 \rangle^{-s} |\mathbf{U}^0|_{H^s}.
$$

The general case $0 \le r \le s$ follows by the interpolation inequality, Proposition A.2, and using once again that
$|U - U_N|_{H^s} \le 4|U^0|_{H^s}$ by the triangle inequality.

Let us now explain why the same result holds in the general case $U^0 \in H^s((2\pi\mathbb{T})^d)^n$. Consider $(U_k^0)_{k \in \mathbb{N}}$ a sequence (constructed by Fourier truncation) such that for all $k \in \mathbb{N}$, $U_k^0 \in H^{s+1}((2\pi\mathbb{T})^d)^n$ and $U_k^0 \to U^0$ in $H^s((2\pi\mathbb{T})^d)^n$ as $k \to \infty$, and $|U_k^0|_{H^s} \le |U^0|_{H^s}$. Then we can apply the above for each $k \in \mathbb{N}$ and infer that $U_k$ (respectively $U_{k,N}$) the solutions to (1.1) (respectively (2.1)) emerging from the initial data $U_k^0$ satisfy for all $0 \le t \le T/|U^0|_{H^s}$ and for any $0 \le r \le s$,

$$
| (U_k - U_{k,N})(t, \cdot) |_{H^r} \le C |U^0|_{H^s} N^{r-s},
$$

where $C > 0$ depends only on $s$ and $M$, and in particular is uniform with respect to $k$. We now pass to the
limit as $k \to \infty$. By standard estimates on the linearized systems from (1.1) and (2.1) (see e.g. [19, Proposition 7.1.8]), we have $U_k \to U$ and $U_{k,N} \to U_N$ in $C([0, T/|U^0|_{H^s}]; L^2((2\pi T)^d)^n)$ as $k \to \infty$, where we denote $U \in C([0, T/|U^0|_{H^s}]; H^s((2\pi T)^d)^n)$ (respectively $U_N \in C([0, T/|U^0|_{H^s}]; H^s((2\pi T)^d)^n)$) the solutions to (1.1) (respectively (2.1)) emerging from the initial data $U^0$, as in the statement of the Proposition. Because the above estimate is uniform with respect to $k$, we infer as desired that the limits satisfy for all $0 \le t \le T/|U^0|_{H^s}$ and for any $0 \le r \le s$,

$$
| (U - U_N)(t, \cdot) |_{H^r} \le C |U^0|_{H^s} N^{r-s}.
$$

Let us now prove the proposition for general $I^*$ compact subset of $I$. Without loss of generality, we can assume
$0 \in I^*$ and we will focus on positive times, $t \in I^* \cap \mathbb{R}^+$, negative times being obtained by time-symmetry. Let
$d/2 + 1 < s_0 < s$. Denote $M^* = 2 \sup_{t \in I^*} |U|_{H^{s_0}}$ and $C^*$ the constant depending on $s_0$, $s$, $M^*$ as in Proposition 2.2
with $s_0$ playing the role of $s$ and $s$ playing the role of $s'$, and $M = 2 |U^0|_{H^s} \exp(C^*M^*T^*)$, where $T^* := \sup(I^*)$.
We set

$$
J_N = \{ t \in I_N \cap \mathbb{R}^+ : \sup_{t' \in [0, t]} |U_N(t', \cdot)|_{H^{s_0}} \le M^* \text{ and } \sup_{t' \in [0, t]} |U_N(t', \cdot)|_{H^s} \le M \},
$$

and our aim is to prove that there exists $N_0 \in \mathbb{N}$ such that for any $N \ge N_0$, $J_N \supset I^* \cap \mathbb{R}^+$. We use the continuity argument. By continuity of $U_N \in C(I_N; H^s((2\pi\mathbb{T})^d)^n)$, we have that $J_N$ is a non-empty closed interval. Let us now prove that $J_N \cap I^* \cap \mathbb{R}^+$ is an open subset of $I^* \cap \mathbb{R}^+$. Let $t \in J_N \cap I^* \cap \mathbb{R}^+$. We can follow the proof of the first part of the proposition replacing the bound $|U|_{H^s} + |U_N|_{H^s} \le 4|U^0|_{H^s}$ with $|U|_{H^s} + |U_N|_{H^s} \le M^* + M$ to infer that there exists $C$, depending only on $s$, $M^*$ and $T^*$ such that for all $0 \le r \le s$,

$$
\sup_{t' \in [0, t]} |U(t', \cdot) - U_N(t', \cdot)|_{H^r} \le C(M^* + M) N^{r-s}.
$$
Applying this estimate with $r = s_0$, it follows that there exists $N_0 \in \mathbb{N}$ depending only on $s_0$, $s$, $M^*$ and $T^*$ such that for any $N \ge N_0$,

$$ \sup_{t' \in [0, t]} |\mathbf{U}_N(t', \cdot)|_{H^{s_0}} \le \sup_{t' \in [0, t]} |\mathbf{U}(t', \cdot)|_{H^{s_0}} + \sup_{t' \in [0, t]} |\mathbf{U}(t', \cdot) - \mathbf{U}_N(t', \cdot)|_{H^{s_0}} \le \frac{2}{3} M^*. $$

Moreover, using the second part of Proposition 2.2 (recall $s_0$ playing the role of $s$ and $s$ playing the role of $s'$), we have

$$ \sup_{t' \in [0, t]} |\mathbf{U}_N(t', \cdot)|_{H^s} \le |\mathbf{U}^0|_{H^s} \exp(C^*M^*t) \le \frac{1}{2}M. $$

This shows, using again the continuity of $\mathbf{U}_N \in C(I_N; H^s((2\pi\mathbb{T})^d)^n)$, that there exists $\delta > 0$ such that $[0, t + \delta] \subset J_N$, and hence that $J_N \cap I^* \cap \mathbb{R}^+$ is a non-empty, connected, closed and open subset of $I^* \cap \mathbb{R}^+$. Hence $J_N \supset I^* \cap \mathbb{R}^+$. Moreover, the desired estimate has been proven along the argument. This concludes the proof. $\square$

## 2.2 Discretization with smooth low-pass filters

Consider now solutions to the systems semi-discretized with a smooth low-pass filter, given by (1.5) and (1.6), and which we recall here for the sake of clarity.

$$ \partial_t \mathbf{U}_N + S_N \left( \sum_{j=1}^{d} A_j(\mathbf{U}_N) \partial_{x_j} \mathbf{U}_N \right) = \mathbf{0}, \quad \mathbf{U}_N|_{t=0} = P_N \mathbf{U}^0, \qquad (2.4) $$

$$ \partial_t \mathbf{U}_N + \sum_{j=1}^{d} (A_j^0 + S_N(A_j^1(\mathbf{U}_N)[\circ])) \partial_{x_j} \mathbf{U}_N = \mathbf{0}, \quad \mathbf{U}_N|_{t=0} = P_N \mathbf{U}^0, \qquad (2.5) $$

where $S_N$ is the smooth low-pass filter described in Section 1, and $A_j(\cdot) = A_j^0 + A_j^1(\cdot)$, $A_j^0 = A_j(\mathbf{0})$.

For symmetric systems, there is no great difference between the semi-discretization with sharp versus smooth low-pass filters. The results are the same, although the proofs must be adapted slightly. We outline the results and proofs, but refer to the previous section for technical details.

Our result regarding existence and boundedness of solutions $\mathbf{U}_N$ to the semi-discrete problems (uniformly with respect to $N$) extends to the case of smooth low-pass filters.

**Proposition 2.5 (Uniform estimates).** *The statement of Proposition 2.2 holds replacing (2.1) with (2.4) or (2.5).*

*Proof.* We consider first the system (2.5). As for the case with the sharp low-pass filter, the key ingredient is an apriori estimate

$$ \frac{d}{dt} |\mathbf{U}_N|_{H^{s'}} \le C(|\mathbf{U}_N|_{H^s}) |\mathbf{U}_N|_{H^s} |\mathbf{U}_N|_{H^{s'}} \qquad (2.6) $$

for $s' \ge s > d/2 + 1$ and $\mathbf{U}_N \in C(I_N; H^s((2\pi\mathbb{T})^d)^n)$ solution to (2.5) on $I_N \subset \mathbb{R}$ open time interval. To show this we first notice that, applying $(\mathrm{Id} - P_N)$ to (2.5) and using that $(\mathrm{Id} - P_N)S_N = 0$, we have $\partial_t(\mathrm{Id} - P_N)\mathbf{U}_N + \sum_{j=1}^d A_j^0 \partial_{x_j}(\mathrm{Id} - P_N)\mathbf{U}_N = 0$ and $(\mathrm{Id} - P_N)\mathbf{U}_N|_{t=0} = 0$. By uniqueness of the solution to this initial-value problem, we infer $\mathbf{U}_N = P_N\mathbf{U}_N$ and hence $\mathbf{U}_N \in C^1(I; H^\sigma((2\pi\mathbb{T})^d)^n)$ for all $\sigma \in \mathbb{R}$. Then we apply $\Lambda^{s'}$ to the system (2.5), denote $\dot{\mathbf{U}}_N := \Lambda^{s'}\mathbf{U}_N$ and infer

$$ \partial_t \dot{\mathbf{U}}_N + \sum_{j=1}^{d} (A_j^0 + S_N^{1/2} A_j^1(\mathbf{U}_N)[S_N^{1/2} \circ]) \partial_{x_j} \dot{\mathbf{U}}_N \\ = -S_N \left( \sum_{j=1}^{d} [\Lambda^{s'}, A_j^1(\mathbf{U}_N)] \partial_{x_j} \mathbf{U}_N \right) - S_N^{1/2} \left( \sum_{j=1}^{d} [S_N^{1/2}, A_j^1(\mathbf{U}_N)] \partial_{x_j} \dot{\mathbf{U}}_N \right). \qquad (2.7) $$
Using that $S_N : L^2 \to L^2$ is bounded, as well as the product, composition and commutator estimates — Propositions A.3, A.4, A.5 and A.6 — the terms on the right-hand side can be estimated in $L^2((2\pi\mathbb{T})^d)^n$ as

$$|\text{RHS}|_{L^2} \le C(|U_N|_{H^s}) |U_N|_{H^s} |\dot{U}_N|_{L^2}.$$

By assumption, $A_j(U) = A_j^0 + A_j(U)$ is a symmetric matrix for all $U \in \mathbb{R}^n$ and $j \in \{1, \dots, d\}$. In particular, this implies that $A_j^0 = A_j(0)$ and hence $A_j^1(U)$ are both symmetric matrices, which in turn implies that the operator $A_j^0 + S_N^{1/2} A_j^1(U_N)S_N^{1/2}$ is symmetric for the $L^2((2\pi\mathbb{T})^d)^n$ inner-product for all $j \in \{1, \dots, d\}$. Arguing in the usual manner using integration by parts, we have

$$\left( \sum_{j=1}^{d} (A_j^0 + S_N^{1/2} A_j^1(U_N)[S_N^{1/2}\circ]) \partial_{x_j} \dot{U}_N, \dot{U}_N \right)_{L^2} = -\frac{1}{2} \sum_{j=1}^{d} \left( [\partial_{x_j}, S_N^{1/2} A_j^1(U_N)S_N^{1/2}] \dot{U}_N, \dot{U}_N \right)_{L^2} \\ \le C(|U_N|_{H^s}) |U_N|_{H^s} |\dot{U}_N|_{L^2}^2,$$

where we used the continuous Sobolev embedding $H^s((2\pi\mathbb{T})^d)^n \subset W^{1,\infty}((2\pi\mathbb{T})^d)^n$ (Proposition A.1) and that $S_N^{1/2} : L^2 \to L^2$ is bounded. Testing the identity (2.7) against $\dot{U}_N$, using the Cauchy–Schwarz inequality on the right-hand side and inserting the two estimates above yields the desired inequality (2.6). With this estimate in hand, the rest of the proof follows exactly the proof of Proposition 2.2.

Considering now the system (2.4), we notice that the inequality (2.6) is obtained as above, using additionally that $S_N A_j^0 = S_N^{1/2} A_j^0 S_N^{1/2}$ is self-adjoint, and the result follows. □

The spectral convergence of solutions $U_N$ to the semi-discrete problems towards the corresponding solution to the continuous problem (1.1) as $N \to \infty$ is the identical in the case of sharp or smooth low-pass filters.

**Proposition 2.6 (Convergence).** *The statement of Proposition 2.4 holds replacing (2.1) with (2.4) or (2.5).*

*Proof.* The proof follows the proof of Proposition 2.4, with one modification. Consider $U_N \in C(I_N; H^s((2\pi\mathbb{T})^d)^n)$ solution to (2.5) (considering instead the solution to (2.4) amounts to replacing $A_j^1(U_N)$ with $A_j(U_N)$ in the right-hand side of the following identity, with no consequence). Then the difference $D_N := U - U_N$ satisfies

$$\partial_t D_N + \sum_{j=1}^{d} A_j(U_N) \partial_{x_j} D_N + \sum_{j=1}^{d} (A_j(U) - A_j(U_N)) \partial_{x_j} U = (\text{Id} - S_N) \left( \sum_{j=1}^{d} A_j^1(U_N) \partial_{x_j} U_N \right).$$

Instead of applying $S_N$ as in the proof of Proposition 2.4, we apply $S_{N/2}$, noting that $S_{N/2}(\text{Id} - S_N) = 0$. Then we have

$$\partial_t S_{N/2} D_N + S_{N/2}^{1/2} \left( \sum_{j=1}^{d} A_j(U_N) \partial_{x_j} S_{N/2}^{1/2} D_N \right) + S_{N/2}^{1/2} \left( \sum_{j=1}^{d} [S_{N/2}^{1/2}, A_j(U_N)] \partial_{x_j} D_N \right) \\ + S_{N/2} \left( \sum_{j=1}^{d} (A_j(U) - A_j(U_N)) \partial_{x_j} U \right) = 0.$$

Because $S_{N/2}$ satisfies the same commutator estimates as $S_N$ and $\|\text{Id} - S_{N/2}^{1/2}\|_{H^s \to L^2} \le \langle N/4 \rangle^{-s}$, we may then proceed exactly as in the proof of Proposition 2.4. □

### 3 Symmetrizable quasilinear systems

In this section, we consider systems that are symmetrizable in the sense of Friedrichs. That is, we assume that there is an open set $U \subset \mathbb{R}^n$ containing the origin and an operator $S(U)$ which is a symmetrizer for the system (1.1):
**Assumption A.2.** There exists $S(\cdot)$ a Friedrichs-symmetrizer for the system (1.1), that is, there exists an open set $U \subset \mathbb{R}^n$ with $0 \in U$ such that for all $U \in \mathcal{U}$, $S(U)$ is real-valued, symmetric positive definite, and for all $j \in \{1, \dots, d\}$, $S(U)A_j(U)$ is symmetric. We assume that all entries of $S(\cdot)$ are polynomial.

**Remark 3.1.** As we assumed in Assumption A.1 that entries of $A_j(U)$ are polynomial for all $j \in \{1, \dots, d\}$, then if Friedrichs-symmetrizers exist it is always possible to select one whose entries are polynomial. Indeed, the assumptions that $S(\cdot)$ is real-valued and symmetric and that for all $j \in \{1, \dots, d\}$, $S(\cdot)A_j(\cdot)$ is symmetric constitutes a system of linear equations for entries of $S$. By considering the corresponding matrix in the field of real rational fractions (since entries of $A_j(U)$ are real polynomials) and performing Gaussian elimination, we see that the system of linear equations can be solved for (non-identically zero) smooth functions if and only if it can be solved for rational fractions, which then can be chosen polynomials after multiplication of all entries by a common multiple of all denominators. Notice however that the domain of hyperbolicity defined as the open set on which $S(\cdot)$ is positive definite depends on the choice of the symmetrizer.

We have the following standard result [2].

**Proposition 3.2 (Well-posedness).** Suppose that the system (1.1) satisfies Assumptions A.1 and A.2. Let $s > 1 + d/2$, $M > 0$ and $K \subset U$ compact. There exists $C > 0$ and $T > 0$ (depending only on $s$, $M$ and $K$) such that for every $U^0 \in H^s((2\pi\mathbb{T})^d)^n$ such that $|U^0|_{H^s} \le M$ and taking values in $K$, there exists a unique $U \in C(I; H^s((2\pi\mathbb{T})^d)^n)$ maximal-in-time classical solution to (1.1) and $U|_{t=0} = U^0$, and moreover the open time interval $I \supset [0, T/|U^0|_{H^s}]$ and for all $t \in [0, T/|U^0|_{H^s}]$,

$$|U|_{H^s} \le C|U^0|_{H^s} \exp(C|U^0|_{H^s}t).$$

## 3.1 Discretization with smooth low-pass filters

We consider the spatial discretization with smooth low-pass filters first, as it is more similar to the previous section. Recall that the spatial discretization with the smooth low-pass filter $S_N$ amounts to solving

$$\partial_t U_N + S_N \left( \sum_{j=1}^{d} A_j(U_N) \partial_{x_j} U_N \right) = 0, \quad U_N|_{t=0} = P_N U^0. \qquad (3.1)$$

As discussed in the introduction, one would typically prefer in practice the variant

$$\partial_t U_N + \sum_{j=1}^{d} (A_j^0 + S_N(A_j^1(U_N)[\circ])) \partial_{x_j} U_N = 0, \quad U_N|_{t=0} = P_N U^0, \qquad (3.2)$$

where $A_j^0 := A_j(0)$ and $A_j^1(\cdot) := A_j(\cdot) - A_j^0$. However, we face a difficulty that while any symmetrizer for the system (1.1), $S(\cdot)$, readily provides a suitable symmetrizer of the semi-discretized system (3.1), such is not the case for system (3.2), and additional assumptions are needed.

**Assumption A.3 (Compatibility of the symmetrizer).** Supposing the Assumptions A.1 and A.2 hold, and decomposing $A_j(\cdot) = A_j^0 + A_j^1(\cdot)$ where $A_j^0 := A_j(0)$ and $S(\cdot) = S^0 + S^1(\cdot)$ where $S^0 := S(0)$, we have

$$\forall j \in \{1, \dots, d\}, \forall U \in \mathcal{U}, \quad S^0 A_j^0, S^0 A_j^1(U) + S^1(U) A_j^0 \text{ and } S^1(U) A_j^1(U) \text{ are symmetric.}$$

**Remark 3.3.** As can be seen by Taylor-expanding $S(U)A_j(U)$ about the origin, Assumption A.3 holds in particular when $A_j(\cdot)$ and $S(\cdot)$ are linear, that is entries of $A_j^1(\cdot)$ and $S^1(\cdot)$ are homogeneous polynomials of degree 1; see also Remark 3.6.

We have the following bound on solutions $U_N$ to the semi-discretized problems (3.1) or (3.2).
**Proposition 3.4** (Uniform estimates). Suppose that the system (1.1) satisfies Assumptions A.1 and A.2. Let $s > d/2 + 1$, $M > 0$ and $K \subset \mathcal{U}$ compact. There exists $C > 0$ and $T > 0$ (depending only on $s$, $M$, and $K$) such that for every $N \in \mathbb{N}$ and for every $\mathbf{U}^0 \in H^s((2\pi\mathbb{T})^d)^n$ such that $|\mathbf{U}^0|_{H^s} \le M$ and taking values in $K$, there exists a unique $\mathbf{U}_N \in C(I_N; H^s((2\pi\mathbb{T})^d)^n)$ maximal-in-time classical solution to (3.1) and $\mathbf{U}|_{t=0} = P_N\mathbf{U}^0$. The open time interval $I_N \supset [0, T/|\mathbf{U}^0|_{H^s}]$ and for all $0 \le t \le T/|\mathbf{U}^0|_{H^s}$,

$$|\mathbf{U}_N|_{H^s} \le C|\mathbf{U}^0|_{H^s} \exp(C|\mathbf{U}^0|_{H^s}t) \le 2C|\mathbf{U}^0|_{H^s}.$$

Moreover, for any $s' \ge s$, one has $\mathbf{U}_N \in C^1(I_N; H^{s'}((2\pi\mathbb{T})^d)^n)$ and for any $0 < T^* \in I_N$, $M^* > 0$ and $K^* \subset \mathcal{U}$ compact such that $\sup_{t \in [0, T^*]} |\mathbf{U}_N(t, \cdot)|_{H^s} \le M^*$ and $\mathbf{U}_N([0, T^*] \times \mathbb{R}^d) \subset K^*$ there exists $C^* > 0$ depending only on $s$, $s'$, $M^*$ and $K^*$ such that for all $0 \le t \le T^*$,

$$|\mathbf{U}_N|_{H^{s'}} \le C^* |\mathbf{U}^0|_{H^{s'}} \exp(C^* M^* t).$$

The same results holds replacing (3.1) with (3.2) if additionally Assumption A.3 holds.

*Proof.* As in Propositions 2.2 and 2.5, the key ingredient is an apriori estimate, and we focus on the derivation of such estimates for solutions $\mathbf{U}_N \in C^1(I_N; H^\sigma((2\pi\mathbb{T})^d)^n)$ for all $\sigma \in \mathbb{R}$ (recall $\mathbf{U}_N = P_N\mathbf{U}_N$ since $S_N = S_N P_N$). We consider system (3.2) which is the most involved, the case of system (3.1) being obtained in the same way. Indeed, while the symmetrizer of the system (1.1), $S(\mathbf{U})$, readily offers a suitable symmetrizer of the semi-discretized system (3.1), for (3.2) we need to consider a modification, namely

$$\tilde{S}(\mathbf{U})[\circ] = S^0 + S_N^{1/2}(S^1(\mathbf{U})S_N^{1/2}[\circ]). \qquad (3.3)$$

Recall $S_N = \text{Diag}(S_N(D))$ and the symbol $S_N$ is nonnegative, hence $S_N^{1/2} = \text{Diag}(S_N(D)^{1/2})$ is well-defined. Let us first prove that for any $K \subset \mathcal{U}$ compact, there exists $0 < \alpha \le \beta < \infty$ such that for all $\mathbf{U} \subset K$ and $\mathbf{V} \in \mathbb{R}^n$ one has

$$\alpha |\mathbf{V}|_{L^2}^2 \le (\tilde{S}(\mathbf{U})\mathbf{V}, \mathbf{V})_{L^2} \le \beta |\mathbf{V}|_{L^2}^2. \qquad (3.4)$$

Note first that the result holds for $S(\mathbf{U})$ by Assumption A.2, and hence for $S^0 = S(0)$ since $0 \in \mathcal{U}$. The result then follows from the identities

$$\begin{aligned} (\tilde{S}(\mathbf{U})\mathbf{V}, \mathbf{V})_{L^2} &= (S^0\mathbf{V}, \mathbf{V})_{L^2} + (S^1(\mathbf{U})S_N^{1/2}\mathbf{V}, S_N^{1/2}\mathbf{V})_{L^2} \\ &= (S^0(\text{Id} - S_N)^{1/2}\mathbf{V}, (\text{Id} - S_N)^{1/2}\mathbf{V})_{L^2} + (S(\mathbf{U})S_N^{1/2}\mathbf{V}, S_N^{1/2}\mathbf{V})_{L^2} \end{aligned}$$

and

$$|\mathbf{V}|_{L^2}^2 = |(\text{Id} - S_N)^{1/2} \mathbf{V}|_{L^2}^2 + |S_N^{1/2} \mathbf{V}|_{L^2}^2.$$

Now we want to show that for any $s' \ge s > 1 + d/2$ and $\mathbf{U}_N$ solution to (3.2) taking values in $K \subset \mathcal{U}$ compact one has

$$\frac{1}{2} \frac{d}{dt} \left( \tilde{S}(\mathbf{U}_N) \Lambda^{s'} \mathbf{U}_N, \Lambda^{s'} \mathbf{U}_N \right)_{L^2} \le C(\mathcal{K}, |\mathbf{U}_N|_{H^s}) |\mathbf{U}_N|_{H^s} |\mathbf{U}_N|_{H^{s'}}^2. \qquad (3.5)$$

Indeed, if (3.5) holds, applying (3.4) in (3.5) and using Grönwall's Lemma yields the following estimate

$$\begin{aligned} \alpha^{1/2}(\mathcal{K}^*) |\mathbf{U}_N(t, \cdot)|_{H^{s'}} &\le \mathcal{F}_{s'}(t) \le \mathcal{F}_{s'}(0) \exp \left( \alpha^{-1}(\mathcal{K}^*) \int_0^t C(\mathcal{K}^*, |\mathbf{U}_N(\tau, \cdot)|_{H^s}) |\mathbf{U}_N(\tau, \cdot)|_{H^s} \, d\tau \right) \\ &\le \beta^{1/2}(\mathcal{K}^*) |\mathbf{U}_N(0, \cdot)|_{H^{s'}} \exp \left( t \times \alpha^{-1}(\mathcal{K}^*) C(\mathcal{K}^*; \mathcal{M}_s(t)) \mathcal{M}_s(t) \right) \end{aligned}$$

where $\mathcal{K}^*$ is compact with $\mathcal{K} \subset \mathcal{U}^* \subset \mathcal{K}^* \subset \mathcal{U}$ and $\mathcal{U}^*$ open, and where we denote

$$\mathcal{F}_{s'}(\tau) := (\tilde{S}(\mathbf{U}_N)\Lambda^{s'}\mathbf{U}_N, \Lambda^{s'}\dot{\mathbf{U}}_N)_{L^2}^{1/2} \Big|_{t=\tau} \quad \text{and} \quad \mathcal{M}_s(t) := \sup_{\tau \in [0,t]} |\mathbf{U}_N(\tau, \cdot)|_{H^s}.$$
The estimate is valid as long as $U_N(\tau, \cdot)$ takes values in $K^*$ for all $\tau \in [0, t]$. This is ensured by

$$|U_N(\tau, \cdot) - U_N^0|_{L^\infty} \le \int_0^t |\partial_t U_N(\tau, \cdot)|_{L^\infty} d\tau \le t \times C(K^*; M_s(t))M_s(t)$$

which follows from the system (3.2) and continuous Sobolev embedding $H^s((2\pi\mathbb{T})^d)^n \subset W^{1,\infty}((2\pi\mathbb{T})^d)^n$ (Proposition A.1). By means of these estimates, we can employ the continuity argument as in the proof of Proposition 2.2 to conclude the proof.

Let us now prove the estimate (3.5). Apply $\Lambda^{s'}$ to the system (3.2) and denote $\dot{U}_N := \Lambda^{s'}U_N$ to infer

$$\partial_t \dot{U}_N + \sum_{j=1}^{d} (A_j^0 + S_N(A_j^1(U_N)[\circ])) \partial_{x_j} \dot{U}_N = -S_N \left( \sum_{j=1}^{d} [\Lambda^{s'}, A_j^1(U_N)] \partial_{x_j} U_N \right).$$

Applying the operator $\tilde{S}(U_N)$ to the system and using the self-adjointness of $\tilde{S}(U_N)$ we get

$$\frac{1}{2} \frac{d}{dt} (\tilde{S}(U_N) \dot{U}_N, \dot{U}_N)_{L^2} = \frac{1}{2} \left( [\partial_t, \tilde{S}(U_N)] \dot{U}_N, \dot{U}_N \right)_{L^2} - \sum_{j=1}^{d} \left( \tilde{S}(U_N)(A_j^0 + S_N(A_j^1(U_N)[\circ])) \partial_{x_j} \dot{U}_N, \dot{U}_N \right)_{L^2} \\ - \sum_{j=1}^{d} \left( \tilde{S}(U_N) S_N([\Lambda^{s'}, A_j^1(U_N)] \partial_{x_j} \dot{U}_N), \dot{U}_N \right)_{L^2}.$$

By using that $\partial_t$ commutes with $S^0$ and $S_N^{1/2}$ we get

$$\left| \left( [\partial_t, \tilde{S}(U_N)] \dot{U}_N, \dot{U}_N \right)_{L^2} \right| = \left| \left( [\partial_t, S^1(U_N)] S_N^{1/2} \dot{U}_N, S_N^{1/2} \dot{U}_N \right)_{L^2} \right| \le C(K, |U_N|_{H^\infty}) |\partial_t U_N|_{L^\infty} |S_N^{1/2} \dot{U}_N|_{L^2}^2.$$

It follows from the system (3.2), continuous Sobolev embedding $H^s((2\pi\mathbb{T})^d)^n \subset W^{1,\infty}((2\pi\mathbb{T})^d)^n$ (Proposition A.1) and product and composition estimates (Propositions A.3 and A.4) that

$$\left| \left( [\partial_t, \tilde{S}(U_N)] \dot{U}_N, \dot{U}_N \right)_{L^2} \right| \le C(K, |U_N|_{H^s}) |U_N|_{H^s} |U_N|_{H^{s'}}^2.$$

For the third term, we use the boundedness of the operators $\tilde{S}(U_N)$ and $S_N$ and the commutator estimate in Proposition A.5 and find

$$\left| \left( \tilde{S}(U_N) S_N \left( [\Lambda^{s'}, A_j^1(U_N)] \partial_{x_j} U_N \right), \dot{U}_N \right)_{L^2} \right| \le C(K, |U_N|_{H^s}) |U_N|_{H^s} |U_N|_{H^{s'}}^2.$$

Estimating the remaining term requires more care, since $\tilde{S}(U_N)$ is not a perfect symmetrizer for the system. Observe

$$\tilde{S}(U_N)(A_j^0 + S_N(A_j^1(U_N))) = S^0 A_j^0 + S_N^{1/2}(S^1(U_N)A_j^0 + S^0 A_j^1(U_N))S_N^{1/2} + S_N S^1(U_N)A_j^1(U_N)S_N \\ + S_N^{1/2} S^0 [S_N^{1/2}, A_j^1(U_N)] - S_N^{1/2}[S_N^{1/2}, S^1(U_N)]S_N A_j^1(U_N) + S_N S^1(U_N)[S_N, A_j^1(U_N)] \\ = \mathrm{Sym}(U_N) + \mathrm{Com}(U_N).$$

Notice that all terms in the first line —the sum being denoted $\mathrm{Sym}(U_N)$— are symmetric operators by Assumption A.3, while the remaining terms —the sum being denoted $\mathrm{Com}(U_N)$— involve commutators between $S_N^{1/2}$ and either $S^1(U_N)$ or $A_j^1(U_N)$. Integration by parts and symmetry considerations implies that

$$\left| \left( \mathrm{Sym}(U_N) \partial_{x_j} \dot{U}_N, \dot{U}_N \right)_{L^2} \right| \le C(K, |U_N|_{H^s}) |U_N|_{H^s} |U_N|_{H^{s'}}^2.$$

For the other terms we use that $S_N^{1/2}$ satisfies the commutator estimate in Proposition A.6 to infer

$$\left| \left( \mathrm{Com}(U_N) \partial_{x_j} \dot{U}_N, \dot{U}_N \right)_{L^2} \right| \le C(K, |U_N|_{H^s}) |U_N|_{H^s} |U_N|_{H^{s'}}^2.$$

Combining all the estimates shows (3.5), which concludes the proof. $\square$

<figure><img src="image_2.png" /></figure>
Having established uniform bounds for solutions $U_N$ to the semi-discretized system (3.1) and (3.2), we now show the convergence towards corresponding solutions of the underlying continuous problem (1.1) as $N \to \infty$.

**Proposition 3.5 (Convergence).** Suppose that the system (1.1) satisfies Assumptions A.1 and A.2. Let $s > d/2 + 1$ and $U^0 \in H^s((2\pi\mathbb{T})^d)^n$ such that $U^0$ takes values into the hyperbolic domain $\mathcal{U}$. Denote $U \in C(I; H^s((2\pi\mathbb{T})^d)^n)$ the maximal-in-time classical solution to (1.1), and $U_N \in C(I_N; H^s((2\pi\mathbb{T})^d)^n)$ the maximal-in-time classical solution to (3.1). For every compact subset $I^* \subset I$, there is an $N_0 \in \mathbb{N}$ and $C^* > 0$, depending only on $s$, $|I^*|$, $\sup_{t \in I^*} |U(t, \cdot)|_{H^s}$ and $K^* \subset \mathcal{U}$ compact such that $U(I^* \times \mathbb{R}^d) \subset K^*$ such that for all $N \ge N_0$, one has $I_N \supset I^*$ and for any $0 \le r \le s$,

$$ \sup_{t \in I^*} |(U - U_N)(t, \cdot)|_{H^r} \le C^* M^* N^{r-s}. $$

The same results holds replacing (3.1) with (3.2) if additionally Assumption A.3 holds.

*Proof.* The proof is similar to that of Propositions 2.4 and 2.6, and we only sketch the main arguments. We consider the system (3.1); system (3.2) can be treated in a similar way after introducing the symmetrizer $\tilde{S}(U)$ defined in (3.3) as in the proof of Proposition 3.4. Denote $D_N := U - U_N$ and use that $(\mathrm{Id} - S_N)S_{N/2}^{1/2} = 0$ to infer

$$ \partial_t S_{N/2}^{1/2} D_N + \sum_{j=1}^{d} A_j(U_N) \partial_{x_j} S_{N/2}^{1/2} D_N + \sum_{j=1}^{d} [S_{N/2}^{1/2}, A_j(U_N)] \partial_{x_j} D_N \\ + S_{N/2}^{1/2} \left( \sum_{j=1}^{d} (A_j(U) - A_j(U_N)) \partial_{x_j} U \right) = 0. $$

Applying the symmetrizer $S(U_N)$ and testing against $S_{N/2}^{1/2} D_N$ yields,

$$ \frac{1}{2} \frac{d}{dt} \left( S(U_N) S_{N/2}^{1/2} D_N, S_{N/2}^{1/2} D_N \right)_{L^2} = \frac{1}{2} \left( [\partial_t, S(U_N)] S_{N/2}^{1/2} D_N, S_{N/2}^{1/2} D_N \right)_{L^2} \\ - \sum_{j=1}^{d} \left( S(U_N) A_j(U_N) \partial_{x_j} S_{N/2}^{1/2} D_N, S_{N/2}^{1/2} D_N \right)_{L^2} \\ - \sum_{j=1}^{d} \left( S(U_N) [S_{N/2}^{1/2}, A_j(U_N)] \partial_{x_j} D_N, S_{N/2}^{1/2} D_N \right)_{L^2} \\ - \sum_{j=1}^{d} \left( S(U_N) S_{N/2}^{1/2} \left( (A_j(U) - A_j(U_N)) \partial_{x_j} U \right), S_{N/2}^{1/2} D_N \right)_{L^2}. $$

The first two terms can be estimated in the standard manner using that $S(\cdot)$ is a symmetrizer, i.e. Assumption A.2. The third term is estimated using the regularizing properties of the commutator with $S_{N/2}^{1/2}$, Proposition A.6. The fourth term is estimated thanks to Assumption A.1. Altogether we obtain

$$ \frac{1}{2} \frac{d}{dt} \left( S(U_N) S_{N/2}^{1/2} D_N, S_{N/2}^{1/2} D_N \right)_{L^2} \le C \left( |U|_{H^s}, |U_N|_{H^s} \right) \left( |U|_{H^s} + |U_N|_{H^s} \right) \left( |S_{N/2}^{1/2} D_N|_{L^2} + |D_N|_{L^2} \right) |S_{N/2}^{1/2} D_N|_{L^2}. $$

With this estimate in hand, we can follow the proof of Proposition 2.4. Using that $\| \mathrm{Id} - S_{N/2}^{1/2} \|_{H^s \to L^2} \le \langle N/4 \rangle^{-s}$, the coercivity of $S(U_N)$ —see (3.4)— and the uniform estimates for $U_N$ stated in Proposition 3.4 we infer (by Grönwall's lemma), denoting $M^* := 2 \sup_{t \in I^*} |U(t, \cdot)|_{H^{s_0}}$ with $d/2 + 1 < s_0 < s$ and assuming $\sup_{t \in I^*} |U_N(t, \cdot)|_{H^{s_0}} \le M^*$ and $U_N(I^* \times \mathbb{R}^n) \subset K^* \subset \mathcal{U}$,

$$ |S_{N/2}^{1/2} D_N(t, \cdot)|_{L^2} \le C(K^*, M^*) M^* \left( |S_{N/2}^{1/2} D_N|_{t=0}|_{L^2} + t N^{-s} |U^0|_{H^s} \right) \exp \left( C(K^*, M^*) M^* t \right), $$
and we have moreover $S_{N/2}^{1/2} \mathbf{D}_N|_{t=0} = 0$ and

$$|\mathbf{D}_N|_{L^2} \le |S_{N/2}^{1/2} \mathbf{D}_N|_{L^2} + \langle N/4 \rangle^{-s} |\mathbf{U}^0|_{H^s}.$$

This yields the desired estimate for $r = 0$, and the general case follows by interpolation. A continuity argument as in Proposition 2.4 allows to secure the bound $\sup_{t \in I^*} |U_N(t, \cdot)|_{H^{s_0}} \le M^*$ and the assumption $U_N(I^* \times \mathbb{R}^n) \subset \mathcal{K}^* \subset \mathcal{U}$ for $N$ sufficiently large (by the convergence $U_N(t, \cdot) \to U(t, \cdot)$ as $N \to \infty$ in $H^{s_0}(\mathbb{R}^d)^n \subset W^{1, \infty}(\mathbb{R}^d)^n$) along the desired estimate, which concludes the proof. $\square$

**Remark 3.6.** Let us comment on the restrictive Assumption A.3 arising in the study of system (3.2). This assumption allows to ensure that we can construct a symmetrizer operator $\tilde{S}(\cdot)$ that is

1. bounded, coercive and self-adjoint for the $L^2(\mathbb{R}^d)^n$ inner-product, and

2. such that $\tilde{S}(U)(A_j^0 + S_N A_j^1(U))$ is self-adjoint up to regularizing operators of order $-1$.

Let us notice that if Assumption A.3 does not hold, it is possible to modify the semi-discretized scheme (3.2) in a way that allows for a symmetrizer operator satisfying at least partially the above requirements. Specifically, consider the system

$$\partial_t U_N + \sum_{j=1}^{d} \sum_{\ell=0}^{m_j} S_N^\ell (A_j^\ell(U_N)[\circ]) \partial_{x_j} U_N = \mathbf{0}, \quad U_N|_{t=0} = P_N U^0, \qquad (3.6)$$

and associated symmetrizer

$$\tilde{S}(U)[\circ] := \sum_{\ell=0}^{m} S_N^{\ell/2} (S^{\ell}(U) S_N^{\ell/2}[\circ]),$$

where we used the convention $S_N^0 = \text{Id}$ and $S_N^\ell = \underbrace{S_N \circ \cdots \circ S_N}_{\ell \text{ times}}$, and decompositions

$$S(U) = \sum_{\ell=0}^{m} S^{\ell}(U), \quad A_j(U) = \sum_{\ell=0}^{m_j} A_j^{\ell}(U) \quad (j \in \{1, \dots, d\})$$

where entries of $S^\ell(\cdot)$ and $A_j^\ell(\cdot)$ are homogeneous polynomials of degree $\ell$.

By Taylor expansion about the origin of $S(\cdot)$ and $S(\cdot)A_j(\cdot)$ and homogeneity we find that for all $U \in \mathcal{U}$, $S^\ell(U)$ and $\sum_{\ell_1+\ell_2=\ell} S^{\ell_1}(U)A_j^{\ell_2}(U)$ are symmetric. This shows that $\tilde{S}(U)$ is self-adjoint for the $L^2(\mathbb{R}^d)^n$ inner-product and, denoting $\tilde{A}_j(U)[\circ] := \sum_{\ell=0}^{m_j} S_N^\ell(A_j^\ell(U)[\circ])$, one has that $\tilde{S}(U)\tilde{A}_j(U)$ is self-adjoint up to regularizing operators of order $-1$. Notice however that some additional restrictions on $U \in \mathcal{U}$ may be necessary to enforce the coercivity of the operator $\tilde{S}(U)$.

## 3.2 Discretization with sharp low-pass filters

In this section we consider the case of spatial discretization through the sharp low-pass filter $P_N = \text{Diag}(P_N(D))$ where $P_N(\cdot) = \mathbf{1}_{[-N,N]^d}(\cdot)$:

$$\partial_t U_N + P_N \left( \sum_{j=1}^{d} A_j(U_N) \partial_{x_j} U_N \right) = \mathbf{0}, \quad U_N|_{t=0} = P_N U^0. \qquad (3.7)$$

The analysis of the previous section fails, due to the lack of good commutator properties of the operator $P_N$. Specifically, when applying a symmetrizer $S(\cdot)$ to the underlying system (1.1) we may write

$$S(U_N) \partial_t U_N + \sum_{j=1}^{d} S(U_N) A_j(U_N) \partial_{x_j} U_N = - \sum_{j=1}^{d} S(U_N) (\text{Id} - P_N) (A_j(U_N) \partial_{x_j} U_N).$$
In order to control the energy functional $\mathcal{F}_s(U_N) := (S(U_N)\Lambda^s U_N, \Lambda^s U_N)^{1/2} \approx |U_N|_{H^s}$ as in Proposition 3.4, we wish to control (uniformly with respect to $N$) the quantity $J_N(U_N, \Lambda^s U_N)$ where

$$J_N(U, V) := (P_N(S(U)(\mathrm{Id} - P_N)(A_j(U)(\partial_{x_j} P_N V))), V)_{L^2}.$$

As we shall see in an example in Section 4, it turns out we cannot improve in general the bound $J_N = O(N|V|_{L^2}^2)$. Notice that in the symmetric cases discussed in Section 2.1, namely when $S = \mathrm{Id}$, we have $J_N = 0$.

In this section we consider symmetrizable systems satisfying the following assumption.

**Assumption A.5.** There exists $S(\cdot)$ and an open set $U \subset \mathbb{R}^n$ with $0 \in U$ such that for all $U \in U$, $S(U)$ is real-valued, symmetric positive definite, and for all $j \in \{1, \dots, d\}$ there exists $S_j^0$ real-valued symmetric matrix with constant coefficients such that for all $U \in U$,

$$A_j(U) = S_j^0 S(U).$$

We assume that all entries of $S(\cdot)$ are polynomial.

**Remark 3.7.** Assumption A.5 is a special case of symmetrizable systems, as it implies Assumptions A.1 and A.2.

**Remark 3.8 (Hamiltonian systems).** Assumption A.5 is motivated by the Hamiltonian structure of the underlying system. Indeed, denote

$$J := \sum_{j=1}^{d} S_j^0 \partial_{x_j}$$

the constant-coefficient skew-symmetric (for the $L^2(\mathbb{R}^d)^n$ inner-product) operator and $\mathcal{H}: U \in U \to \mathbb{R}$ coercive functional such that for all $U \in U$,

$$\mathrm{Hess}(\mathcal{H}(U)) = S(U).$$

Then we remark that under the Assumption A.5 (1.1) takes the Hamiltonian form

$$\partial_t U + J(\nabla_U \mathcal{H}(U)) = 0,$$

where $\nabla_U \mathcal{H}: U \to \mathbb{R}^n$ is the Jacobian of $\mathcal{H}$.

Moreover, noticing that (3.7) also enjoys a Hamiltonian structure,

$$\partial_t U_N + J(\nabla_U \mathcal{H}_N(U_N)) = 0$$

where $\mathcal{H}_N(U) = \mathcal{H}(P_N U)$, we find that $\mathrm{Hess}(\mathcal{H}_N(U_N)) = P_N S(U_N) P_N$ is a symmetrizer of the system (3.7).

Under this assumption, we have the following bound on solutions $U_N$ to the semi-discretized problems (3.7).

**Proposition 3.9 (Uniform estimates).** Under the Assumption A.5, the statement of Proposition 3.4 holds replacing (3.1) with (3.7).

*Proof.* We follow very closely the proof of Proposition 3.4, and only sketch how the necessary estimates can be obtained. Apply $P_N S(U_N) P_N \Lambda^{s'}$ to the system (3.7) and use the identity $A_j(U_N) = S_j^0 S(U_N)$ to infer

$$P_N(S(U_N)P_N(\partial_t \Lambda^{s'} U_N)) + \sum_{j=1}^{d} P_N(S(U_N)(P_N S_j^0 P_N(S(U_N)\partial_{x_j} \Lambda^{s'} U_N))) \\ = -P_N\left(\sum_{j=1}^{d} S(U_N)P_N([\Lambda^{s'}, S_j^0 S(U_N)]\partial_{x_j} U_N)\right)$$
where we used that $P_N^2 = P_N$ commutes with $S_j^0$ and $\Lambda^{s'}$. We can now test the identity against $\Lambda^{s'}U_N$ and use the self-adjointness of $P_N$, $S_j^0$ and $S(U_N)$ as well as the identity $U_N = P_N U_N$ to infer

$$ \frac{1}{2} \frac{d}{dt} (S(U_N)\Lambda^{s'} U_N, \Lambda^{s'} U_N)_{L^2} = \frac{1}{2} ([\partial_t, S(U_N)]\Lambda^{s'} U_N, \Lambda^{s'} U_N)_{L^2} \\ + \sum_{j=1}^{d} (S_j^0 P_N ([\partial_{x_j}, S(U_N)]\Lambda^{s'} U_N), P_N (S(U_N)\Lambda^{s'} U_N))_{L^2} \\ - \sum_{j=1}^{d} (S(U_N)P_N ([\Lambda^{s'}, S_j^0 S(U_N)]\partial_{x_j} U_N), \Lambda^{s'} U_N)_{L^2}. $$

We then proceed as in Proposition 3.4 and obtain the energy estimate valid as long as $U_N$ takes values into $\mathcal{K} \subset \mathcal{U}$ compact:

$$ \frac{1}{2} \frac{d}{dt} (S(U_N)\Lambda^{s'} U_N, \Lambda^{s'} U_N)_{L^2} \le C |U_N|_{H^s} |U_N|_{H^{s'}}^2, $$

where the constant $C$ depends only on $s, s', \mathcal{K}$ and non-decreasingly on $|U_N|_{H^s}$. We also have immediately the coercivity of $S(U)$: for any $U \subset \mathcal{K}$ and $V \in \mathbb{R}^n$ one has

$$ \alpha |V|_{L^2}^2 \le (S(U)V, V) \le \beta |V|_{L^2}^2, $$

where $0 < \alpha \le \beta < \infty$ depend uniquely on $\mathcal{K}$. These two ingredients yield the desired result. $\square$

Having established uniform bounds for solutions $U_N$ to the semi-discretized system (3.7), we infer the convergence towards corresponding solutions of the underlying continuous problem (1.1) as $N \to \infty$.

**Proposition 3.10 (Convergence).** Under Assumption A.5, the statement of Proposition 3.5 holds replacing (3.1) with (3.7).

*Proof.* The proof is identical to that of Proposition 3.5. $\square$

# 4 Numerical experiments for the Saint-Venant system

We shall illustrate our findings and investigate numerically the standard Saint-Venant (or shallow water) system

$$ \left\{ \begin{array}{l} \partial_t \eta + \nabla \cdot ((1 + \eta)\mathbf{u}) = 0, \\ \partial_t \mathbf{u} + \nabla \eta + (\mathbf{u} \cdot \nabla)\mathbf{u} = \mathbf{0}, \end{array} \right. \qquad (4.1) $$

which describes the propagation of shallow water waves in the flat-bottom situation; see [17]. Specifically, the scalar variable $\eta$ describes the elevation of the surface of a layer of homogeneous, incompressible and inviscid fluid and the variable $\mathbf{u}$ represents the layer-averaged horizontal velocity of fluid particles (both depending on time and horizontal space). The gravitational constant and reference depth have been set to $g = 1$ and $H = 1$.

It will be interesting to consider the following variant (when $d = 2$)

$$ \left\{ \begin{array}{l} \partial_t \eta + \nabla \cdot ((1 + \eta)\mathbf{u}) = 0, \\ \partial_t \mathbf{u} + \nabla \eta + \frac{1}{2} \nabla (|\mathbf{u}|^2) = \mathbf{0}. \end{array} \right. \qquad (4.2) $$

While the two systems are identical when $d = 1$, only the second has a Hamiltonian structure when $d = 2$; see below. As exhibited in the following section, the hyperbolicity domain of the Hamiltonian system (4.2) is a strict subset of the hyperbolicity domain of the standard system (4.1).

Thanks to these features, numerical experiments on the Saint-Venant systems allow to showcase our numerical findings, that we summarize here for the sake of readability.
* Numerical experiments validate our results concerning the spectral convergence of the (semi-)discretized solutions as $N \to \infty$, when Assumptions A.1, A.2 and A.3 hold and smooth low-pass filters are used.

* In the case of sharp low-pass filters, we have not been able to observe numerical instabilities when Assumptions A.1, A.2 and A.3 hold but Assumption A.5 fails.

* Contrarily to sharp low-pass filters, smooth low-pass filters are able to instate a form of stability even outside the domain of hyperbolicity, that is when Assumption A.2 fails.

## 4.1 Analysis of the Saint-Venant system

We can apply the analysis of the previous section to systems (4.1) and (4.2) due to the following result.

**Lemma 4.1.** System (4.1) is a symmetrizable hyperbolic system in the sense of Assumptions A.1 and A.2 with hyperbolic domain $\mathcal{U} := \{(\eta, \mathbf{u}) \in \mathbb{R}^{1+d} : 1 + \eta > 0\}$ and

$$S((\eta, \mathbf{u})) = \begin{pmatrix} 1 & \mathbf{0}^{\top} \\ \mathbf{0} & (1 + \eta) \mathrm{Id} \end{pmatrix} \qquad (4.3)$$

where $\mathrm{Id}$ is the identity matrix in $\mathbb{R}^d$. Moreover, the additional Assumption A.3 holds (see Remark 3.3).

System (4.2) satisfies Assumption A.5 (and hence Assumptions A.1 and A.2; see Remark 3.7) with hyperbolic domain $\mathcal{U}_\mathcal{H} := \{(\eta, \mathbf{u}) \in \mathbb{R}^{1+d} : 1 + \eta - |\mathbf{u}|^2 > 0\}$ and

$$S_{\mathcal{H}}((\eta, \mathbf{u})) = \begin{pmatrix} 1 & \mathbf{u}^{\top} \\ \mathbf{u} & (1 + \eta) \mathrm{Id} \end{pmatrix} \quad \text{and} \quad S_j^0 = \begin{pmatrix} 0 & \mathbf{e}_j^{\top} \\ \mathbf{e}_j & \mathbf{0} \end{pmatrix} \quad \text{where } \mathbf{e}_1 = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \mathbf{e}_2 = \begin{pmatrix} 0 \\ 1 \end{pmatrix} \qquad (4.4)$$

(when $d = 2$, set $\mathbf{e}_j = 1$ for the analogous definitions for $d = 1$) and Hamiltonian energy

$$\mathcal{H}((\eta, \mathbf{u})) = \frac{1}{2} \int_{(2\pi\mathbb{T})^d} \eta^2 + (1 + \eta) |\mathbf{u}|^2 d\mathbf{x}.$$

**Remark 4.2.** As aforementioned, systems (4.1) and (4.2) are identical when $d = 1$, and hence enjoy both properties. Notice that, when $d = 2$, the domain of hyperbolicity of system (4.2), $\mathcal{U}_\mathcal{H}$, is strictly embedded in the domain of hyperbolicity of system (4.1), $\mathcal{U}$, while only the former satisfies Assumption A.5, associated with its Hamiltonian formulation.

*Proof.* The systems (4.1) and (4.2) can be reformulated as

$$\partial_t \mathcal{U} + \sum_{j=1}^{d} A_j(\mathcal{U}) \partial_{x_j} \mathcal{U} = \mathbf{0}, \qquad (4.5)$$

with $\mathcal{U} = (\eta, \mathbf{u})$ and

$$A_j((\eta, \mathbf{u})) = \begin{pmatrix} u_j & (1 + \eta)\mathbf{e}_j^{\top} \\ \mathbf{e}_j & u_j \mathrm{Id} \end{pmatrix}$$

for system (4.1) and

$$A_j((\eta, \mathbf{u})) = \begin{pmatrix} u_j & (1 + \eta)\mathbf{e}_j^{\top} \\ \mathbf{e}_j & \mathbf{e}_j \mathbf{u}^{\top} \end{pmatrix}$$

for system (4.2). It is then straightforward to check the assumptions. $\square$

<figure><img src="image_3.png" /></figure>
Recall the spatial discretization of the system (4.5) with the smooth low-pass filter $S_N$,

$$ \partial_t U_N + S_N \left( \sum_{j=1}^{d} A_j(U_N) \partial_{x_j} U_N \right) = 0, \quad U_N|_{t=0} = P_N U^0, \qquad (4.6) $$

or

$$ \partial_t U_N + \sum_{j=1}^{d} (A_j^0 + S_N(A_j^1(U_N)[\circ])) \partial_{x_j} U_N = 0, \quad U_N|_{t=0} = P_N U^0, \qquad (4.7) $$

and the spatial discretization with the sharp low-pass filter $P_N$,

$$ \partial_t U_N + P_N \left( \sum_{j=1}^{d} A_j(U_N) \partial_{x_j} U_N \right) = 0, \quad U_N|_{t=0} = P_N U^0. \qquad (4.8) $$

It follows immediately from Lemma 4.1 that we have convergence of the numerical scheme in all three cases.

**Proposition 4.3. (Convergence)** The statement of Proposition 3.5, concerning spectral convergence of solutions to the semi-discrete systems (4.6) and (4.7), holds for the Saint-Venant system (4.1) whenever $U^0$ takes values in $U := \{(\eta, \mathbf{u}) \in \mathbb{R}^{1+d} : 1 + \eta > 0\}$. If additionally $U^0$ takes values in $U_{\mathcal{H}} := \{(\eta, \mathbf{u}) \in \mathbb{R}^{1+d} : 1 + \eta - |\mathbf{u}|^2 > 0\}$, then Proposition 3.5 also holds for the Hamiltonian Saint-Venant system (4.2).

Whenever $U^0$ takes values in $U_{\mathcal{H}}$, we furthermore have that the statement of Proposition 3.10, concerning spectral convergence of solutions to the semi-discrete system (4.8), holds for the system (4.2).

*Proof.* Lemma 4.1 ensures that the systems (4.1) and (4.2) satisfy the assumptions of Proposition 3.5 with symmetrizers $S(U)$ and $S_{\mathcal{H}}(U)$ respectively, and that (4.2) additionally satisfies the assumptions of Proposition 3.10. The domains $U$ and $U_{\mathcal{H}}$ correspond respectively to the domains for which the symmetrizers $S(U)$ and $S_{\mathcal{H}}(U)$ in (4.3), (4.4) are positive definite, as required in Assumptions **A.2** and **A.5**. $\square$

**Remark 4.4.** As discussed in the previous section, Section 3.2, we require more stringent structural assumptions to show convergence for symmetrizable systems when discretizing with the sharp low-pass filter $P_N$. The Saint-Venant system when $d = 1$ illustrates that we may also have to impose more stringent restrictions on the initial data (namely $U^0$ taking values in $U_{\mathcal{H}}$), even for systems that satisfy the structural assumptions.

This is because any symmetrizer for the underlying system in the sense of Assumption **A.2** satisfying the additional compatibility Assumption **A.3** can be used to construct a symmetrizer to the semi-discrete systems (4.6) and (4.7). For the semi-discretization with $P_N$ (4.8) on the other hand, we use the symmetrizer directly related to the structure of the system through Assumption **A.5**.

Let us illustrate the discussion in Remark 4.4. As discussed in the beginning of Section 3.2, considering a semi-discrete system with sharp low-pass filter emanating from a symmetrizable continuous system with symmetrizer $S(U)$, one wishes to control the energy functional $\mathcal{F}_s(U_N) = (S(U_N)\Lambda^s U_N, \Lambda^s U_N)^{1/2} \approx |U_N|_{H^s}$, which in turn requires to control (uniformly with respect to $N$) the quantity $J_N(U_N, \Lambda^s U_N)$ where

$$ J_N(U, V) := (P_N(S(U)(\mathrm{Id} - P_N)(A_j(U)(\partial_{x_j} P_N V))), V)_{L^2}. $$

In the specific case of the Saint-Venant system (4.1) when $d = 1$, one has

$$ \mathcal{F}_s(U) = \int_{(2\pi\mathbb{T})} (\Lambda^s \eta)^2 + (1 + \eta)(\Lambda^s u)^2 dx $$

and

$$ S(U) = \begin{pmatrix} 1 & 0 \\ 0 & 1 + \eta \end{pmatrix}, \quad A(U) = \begin{pmatrix} u & 1 + \eta \\ 1 & u \end{pmatrix}, \quad U = \begin{pmatrix} \eta \\ u \end{pmatrix}. $$
Let $U := (\eta_p, u_p)$ where $\eta_p(x) := -\frac{1}{2}\cos(px)$, $u_p(x) := \sin(px)$ and $V_N := (0, v_{N,q})$ where $v_{N,q}(x) := \sin((N-q)x)$ with $0 \le q < p \ll N$. A direct calculation yields

$$(\mathrm{Id} - \mathrm{P}_N)(A(U)(\partial_x \mathrm{P}_N \mathbf{V}_N)) = \frac{N-q}{4} \begin{pmatrix} -\cos((N-q+p)x) \\ 2\sin((N-q+p)x) \end{pmatrix},$$

so that

$$\mathrm{P}_N(S(U)(\mathrm{Id} - \mathrm{P}_N)(A(U)(\partial_{x_j} \mathrm{P}_N \mathbf{V}_N))) = \frac{N-q}{8} \begin{pmatrix} 0 \\ -\sin((N-q)x) \end{pmatrix},$$

and hence

$$J_N(U, V_N) = -\frac{\pi}{8}(N-q).$$

This shows that one cannot propagate for positive time (at least in a direct manner) a uniform-in-$N$ control of the energy functional $\mathcal{F}_s(U_N)$ for $U_N$ the solution emerging from initial data $U_N^0 := U + V_N/(N-q)^s$, despite the fact that $U_N^0 \in \mathcal{U}$ since $1+\eta_N^0 \ge 1/2 > 0$, and $|U_N^0|_{H^s} \approx 1$. Notice also that $U_N^0 \notin \mathcal{U}_\mathcal{H}$ since $1+\eta_N^0(\frac{\pi}{2p})-|u_N^0(\frac{\pi}{2p})|^2 = 0$, but one could enforce $U_N^0 \in \mathcal{U}_\mathcal{H}$ while keeping valid all previous statements by considering e.g. $u_p(x) = \frac{1}{2}\sin(px)$. In that case, the Hamiltonian structure allows to propagate the functional $\mathcal{F}_{\mathcal{H},s}(U) \approx |U_N|_{H^s}$ with

$$\mathcal{F}_{\mathcal{H},s}(U) := (S_{\mathcal{H}}(U)\Lambda^s U, \Lambda^s U)^{1/2} = \int_{(2\pi\mathbb{T})} (\Lambda^s \eta)^2 + (1+\eta)(\Lambda^s u)^2 + 2u(\Lambda^s \eta)(\Lambda^s \eta) \,dx.$$

## 4.2 Numerical experiments in dimension one

We seek numerical approximations to (4.1) with $d=1$ (or, equivalently, (4.2)), $\eta_M, u_M$, in terms of finite Fourier sums of the form

$$f(x) = \sum_{k=-M+1}^{M} a_k \exp(ikx),$$

and similarly for $u_M$. The vectors $\eta = (\eta_M(x_1), \dots, \eta_M(x_{2M}))$, $u = (u_M(x_1), \dots, u_M(x_{2M}))$ contain the values of $\eta_M, u_M$ at regularly spaced collocation points $x_n = -\pi + \pi n/M, n = 1, \dots, 2M$. We use the discrete Fourier Transform, computed efficiently with a Fast Fourier transform (FFT), to find

$$\eta_M(x) = \sum_{k=-M+1}^{M} \hat{\eta}_k \exp(ikx), \quad u_M(x) = \sum_{k=-M+1}^{M} \hat{u}_k \exp(ikx),$$

where $\hat{\eta} = (\hat{\eta}_{-M+1}, \dots, \hat{\eta}_M)$, $\hat{u} = (\hat{u}_{-M+1}, \dots, \hat{u}_M)$ are the coefficients of the Fast Fourier transform of $\eta, u$, and $\hat{k} = (k/\pi : k = -M+1, \dots, M)$ are the discrete Fourier modes. Abusing notation, we will incorrectly refer to $\hat{\eta}_k, \hat{u}_k$ as Fourier coefficients (they are related to the coefficients $c_j$ of infinite Fourier series $f(x) = \sum_{j \in \mathbb{Z}} c_j \exp(i jx)$ through $\hat{f}_k = \sum_{j \in \mathbb{Z}} c_{k+2jM})$. For functions $f \in H^s(2\pi\mathbb{T})$, the error due to this aliasing effect is of order $O(M^{-s})$.

Spatial differentiation is now obtained by multiplying the Fourier coefficients with $ik$. Nonlinear operations are computed pointwise on collocation points $x_n$, via inverse Fast Fourier transform. This procedure leads, in general, to aliasing errors. For polynomial nonlinearities (such as for the Saint-Venant system), one can use so-called dealiasing techniques to remove these errors. For quadratic nonlinearities, one may for example use Orszag's 3/2-rule [20], which consists in adding a sufficient number of Fourier modes with coefficients set to zero. For more information on spectral methods and dealiasing techniques, we refer to [5] and [25].

In our numerical codes, to remove aliasing errors from the nonlinear terms while still working with vectors $\hat{\eta}, \hat{u}$ of fixed length, we shall set the highest 1/3 of the Fourier modes to zero. As we numerically compute approximate solutions of the semi-discretized equations (4.7) and (4.8), this procedure is naturally performed when applying sharp or smooth low-pass filters, $P_N, S_N$, with $N < 2M/3$. For the smooth low-pass filter, we use the example from the introduction, that is $S_N = \mathrm{Diag}(S_N(D))$ with $S_N(\cdot) = S(\cdot/N)$ and $S(\cdot) = \max(0, \min(1, 2-2|\cdot|))^2$. For the
<figure><img src="image_4.png" /><figcaption>(a) Plot of the initial data (4.9), where $\eta^0$ is the initial surface profile, and $u^0$ is the initial velocity.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Decay of the Fourier coefficients for the initial data (4.9). The blue and orange points show $E_s(U_N)|_{t=0}$ for $s = 0, 1$ respectively. To illustrate, the blue and orange lines have slopes $-2$ and $-1$ respectively.</figcaption></figure>

<caption>Figure 1: Experiments with initial data (4.9).</caption>

discretization with the smooth low-pass filter, we shall consider only the version (4.7) where the low-pass filter is only applied to nonlinear terms.

This procedure of semi-discretization in space yields a system of differential equations in time for the Fourier coefficients $\hat{\eta}$, $\hat{u}$. We approximately solve this initial-value problem using an explicit Runge-Kutta 4 method. All numerical simulations are made using the Julia package WaterWaves1D [9] and can be reproduced using the scripts available at WaterWaves1D.jl/examples/StudySaintVenant.jl.

From now on, we denote the number of collocation points by $2M$, and let $N = [2M/3]$, that is, the greatest integer smaller than $2M/3$. Our numerical scheme maintains the highest $1/3$ of the Fourier modes to zero at each time-step, see the discussion above. Abusing notation, we will refer to the fully-discretized numerical solution as $U_N = (\eta_N, u_N)$, since only $2N$ Fourier modes are nonzero. This convention means that $N$ plays the same role in this section as in the previous, analytical sections. We will compute the numerical solution with $2M = 2^j, j = 6, \dots, 15$ collocation points and use time step $dt = 10^{-5}$. The time step is an order of magnitude smaller than needed to avoid stability issues, and small enough to ensure the error due to the spatial discretization dominates. We will use the solution computed with $2M = 2^{15}$ and sharp low-pass filter as a reference solution $U_{\text{ref}} = (\eta_{\text{ref}}, u_{\text{ref}})$, and compute the relative error of the numerical solutions $U_N = (\eta_N, u_N)$ by comparing with the reference solution:

$$E_s(U_N) = \frac{|U_N - U_{\text{ref}}|_{H^s}}{|U_{\text{ref}}|_{H^s}}.$$

The norms will be computed approximately using the Fourier coefficients of the numerical solutions.

We solve numerically the Saint-Venant system (4.1) in one spatial dimension. For the tested initial data in $\mathcal{U}_{\mathcal{H}}$, numerical results are in agreement with the analysis. To study the experimental convergence, we consider the following initial data for $\alpha > 0$,

$$\eta^0(x) = \frac{1}{2} \exp(-|x|^\alpha) \exp(-4x^2), \quad u^0(x) = 0. \qquad (4.9)$$

Notice $U^0 = (\eta^0, u^0)$ satisfies both $1 + \eta^0 > 0$ and the stricter condition $1 + \eta^0 - (u^0)^2 > 0$. The initial surface is a heap of water situated at the origin. Both $\eta^0, u^0$ decay to machine precision near $-\pi, \pi$ and can therefore be seen as periodic. Moreover, $U^0 \in H^{\alpha+1/2}(2\pi\mathbb{T})^2$. We let $\alpha = 1.5$ and simulate the time-evolution up to a final time $T = 0.5$ with either sharp or smooth lowpass filters applied to the nonlinear terms. A plot of the initial data as well as the decay of its Fourier coefficients (through $E^s(U_N)|_{t=0} = \frac{|(\text{Id} - P_N)U_{\text{ref}}|_{t=0}|_{H^s}}{|U_{\text{ref}}|_{t=0}|_{H^s}}$) is shown in Figure 1.

Figure 2a show log-log plots of the error $E_s$ at time $T = 0.5$ for the numerical solution with $2M = 2^j$ where
<figure><img src="image_6.png" /></figure>

(a) Plot illustrating the convergence of the numerical schemes (4.7) and (4.8) as the number of collocation points $2M$ increases. The plot shows the relative error of the numerical solution for initial data (4.9) in $H^2(2\pi\mathbb{T})^2$ measured in the $L^2$-norm, $E_0$ and in the $H^1$-norm, $E_1$ for $2M = 2^j, j = 6, \dots, 14$ when using either sharp or smooth low-pass filters. To illustrate, the blue and orange lines have slopes $-2$ and $-1$ respectively. The numerical scheme exhibits spectral convergence with both sharp and smooth low-pass filters.

<table><thead><tr><th rowspan="2">2M</th><th colspan="2">Sharp low-pass filter</th><th colspan="2">Smooth low-pass filter</th></tr><tr><th>EOC<sub>0</sub></th><th>EOC<sub>1</sub></th><th>EOC<sub>0</sub></th><th>EOC<sub>1</sub></th></tr></thead><tbody><tr><td>2<sup>6</sup></td><td>1.69</td><td>0.7</td><td>2.75</td><td>1.7</td></tr><tr><td>2<sup>7</sup></td><td>2.02</td><td>1.01</td><td>1.93</td><td>1.00</td></tr><tr><td>2<sup>8</sup></td><td>2.02</td><td>1.01</td><td>1.97</td><td>0.99</td></tr><tr><td>2<sup>9</sup></td><td>2.04</td><td>1.03</td><td>1.95</td><td>0.97</td></tr><tr><td>2<sup>10</sup></td><td>2.03</td><td>1.02</td><td>1.94</td><td>0.96</td></tr><tr><td>2<sup>11</sup></td><td>2.03</td><td>1.02</td><td>1.98</td><td>0.98</td></tr><tr><td>2<sup>12</sup></td><td>2.01</td><td>0.99</td><td>1.98</td><td>0.99</td></tr><tr><td>2<sup>13</sup></td><td>1.95</td><td>0.92</td><td>1.98</td><td>0.98</td></tr></tbody></table>

(b) Experimental order of convergence for the numerical solution with initial data (4.9) for both sharp and smooth low-pass filters. The Experimental order of convergence is measured in the $L^2$-norm, $EOC_0$, and in the $H^1$-norm, $EOC_1$.

<caption>Figure 2: Experiments with initial data (4.9).</caption>

$j = 6, \dots, 14$ computed using sharp and smooth low-pass filters. Figure 2a shows the relative error measured in the $L^2$-norm, $E_0(U_N)$ and in the $H^1$-norm, $E_1(U_N)$.

The experimental order of convergence, given by

$$ EOC_s = \frac{\log(E_s(U_N)/E_s(U_{2N}))}{\log(2)}, $$

is given in Figure 2b. With the initial data in $H^2(2\pi\mathbb{T})^d$, we expect from our analysis that the $L^2$-error should decay as $O(N^{-2})$ and the $H^1$-error should decay as $O(N^{-1})$ when using both the sharp and smooth low-pass filter. This aligns with our numerical results. The absolute error is slightly larger when using the smooth low-pass filter, which is to be expected since applying $S_N$ removes more information at each time step than does $P_N$. Taking as initial data (4.9) with other values of $\alpha$ (we have tested $\alpha = 1, 2.5, 3$) also yields the expected results.

Let us now consider initial data $(\eta^0, u^0) \in \mathcal{U} \setminus \mathcal{U}_\mathcal{H}$, that is satisfying $1 + \eta^0 > 0$, but not the stricter condition $1 + \eta^0 - (u^0)^2 > 0$. In particular, we will consider the initial data from the example at the end of Section 4.1, with $q = 0, p = 1, s = 2$:

$$ U_N^0 = (\eta^0, u_N^0), \quad \eta_N^0(x) = -\frac{1}{2} \cos(x), \quad u_N^0(x) = \sin(x) + \frac{\sin(Nx)}{N^2} \quad (4.10) $$

Notice $f(x) = 1 - \frac{1}{2}\cos(x) - \sin^2(x) < 0$ for $x \in (-\frac{\pi}{2}, -\frac{\pi}{3}) \cup (\frac{\pi}{3}, \frac{\pi}{2})$. On the other hand, $1 + \eta^0(x) > 0$ for all $x \in 2\pi\mathbb{T}$. We therefore expect the numerical scheme to converge when using the smooth low-pass filter, but that instabilities may emerge when discretizing with the sharp low-pass filter. Convergence plots for the numerical solutions in both cases are shown in Figure 3b. In line with the analysis, we have convergence when using the smooth low-pass filter. However, this is also true for the sharp low-pass filter. Despite many attempts, we have not been able to observe numerical instabilities for initial data we have tested satisfying $1 + \eta^0 > 0$ but violating $1 + \eta^0 - (u^0)^2 > 0$.

Interestingly, we do observe a difference between sharp and smooth low-pass filters when (barely) violating the
<figure><img src="image_7.png" /><figcaption>(a) Plot of the initial data (4.10), where $\eta^0$ is the initial surface profile, and $u^0$ is the initial velocity.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Plot illustrating the convergence of the numerical schemes (4.7) and (4.8) as the number of collocation points $2M$ increases. The plot shows the relative error of the numerical solution measured in the $L^2$-norm, $E_0$ and in the $H^1$-norm, $E_1$ for $2M = 2^j, j = 6, \dots, 14$ when using either sharp or smooth low-pass filters. To illustrate, the blue and orange lines have slopes $-2$ and $-1$ respectively. The numerical scheme exhibits spectral convergence with both sharp and smooth low-pass filters.</figcaption></figure>

<caption>Figure 3: Experiments with initial data (4.10).</caption>

non-cavitation assumption $1 + \eta^0 > 0$. Let

$$ \eta_{0,N}(x) = -\cos(x), \quad u_{0,N}(x) = \sin(x) + \frac{\sin(Nx)}{N^2}. \qquad (4.11) $$

Figure 4 shows the second derivative of the velocity at time $T = 0.1$ for smooth and sharp low-pass filters for $2M = 2^{10}$ and $2M = 2^{12}$. The second derivative is uniformly bounded when using the smooth low-pass filter whereas it is not around the point $x = 0$ where the non-cavitation assumption is violated when using the sharp low-pass filter.

## 4.3 Numerical experiments in dimension two

The numerical simulations of the Saint-Venant system with $d = 2$ are analogous to the case when $d = 1$ and have been executed using the same Julia package WaterWaves1D [9]. They are also reproducible using the scripts available at WaterWaves1D.jl/examples/StudySaintVenant.jl. We let $2M$ denote the number of collocation points in each of the two spatial dimensions, $x$ and $y$, which form a grid with $4M^2$ collocation points. As in the previous section, we set $N = [2M/3]$ and let $P_N, S_N$ be as described in the introduction. In particular, notice $S_N$ is now a composition of one-dimensional low-pass filters: $S_N := \text{Diag}(S_N(D))$ with $S_N((k_1, k_2)) := S(k_1/N)S(k_2/N)$, where we set $S(\cdot) := \max(0, \min(1, 2-2|\cdot|))^2$. We denote numerical approximations by $\eta_N, u_N, v_N$. The values at collocation points $\eta, u, v$ and associated discrete Fourier modes $\hat{\eta}, \hat{u}, \hat{v}$ are $2M$-by-$2M$ matrices.

Our main interest in studying numerically the systems in two dimensions is to examine whether we observe any difference between the Hamiltonian and non-Hamiltonian version of the Saint-Venant system, respectively (4.2) and (4.1). While there is indeed a difference with respect to stability of the numerical schemes—we observe instabilities with the sharp low-pass filter when violating $1 + \eta^0 - (u^0)^2 - (v^0)^2 > 0$ for (4.2), but not for (4.1)— this can be explained by the difference between the hyperbolicity domains $\mathcal{U}_\mathcal{H}$ and $\mathcal{U}$ rather than by the presence or absence of a Hamiltonian structure. Because the setting of dimension $d = 2$ is computationally costlier, we use timestep $dt = 5 \times 10^{-4}$ in our numerical experiments and, when calculating relative errors, we take as a reference solution the numerical solution computed with $2M = 2^{10}$.
<figure><img src="image_9.png" /><figcaption>(a) Plot of the second derivative of the velocity $u_N$ at time $T = 0.1$ with initial data (4.11) with $2M = 2^{10}$. The orange and blue lines show the solutions found by using the numerical schemes (4.7) and (4.8) with smooth and sharp low-pass filters respectively.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Plot of the second derivative of the velocity $u_N$ at time $T = 0.1$ with initial data (4.11) with $2M = 2^{12}$. The orange and blue lines show the solutions found by using the numerical schemes (4.7) and (4.8) with smooth and sharp low-pass filters respectively.</figcaption></figure>

<caption>Figure 4: Experiments with initial data (4.11).</caption>

We test the numerical method on initial data in $H^s((2\pi\mathbb{T})^2)^3$ of the form

$$ \begin{aligned} \eta^0(x, y) &= (h_0 - 1) \cos(x) \cos(y), \\ u^0(x, y) &= u_1 \sin(x) \cos(y) + u_h \frac{\sin(Nx) \cos(Ny)}{N^s}, \\ v^0(x, y) &= v_1 \cos(x) \sin(y) + v_h \frac{\cos(Nx) \sin(Ny)}{N^s}, \end{aligned} \quad (4.12) $$

where $h_0 > 0, u_1, v_1, s \ge 0$ are real numbers. When $u_1 = v_1$ and $u_h = v_h$, the initial data is irrotational, and systems (4.1) and (4.2) are equivalent.

For the standard, non-Hamiltonian Saint-Venant system (4.1) our numerical results when $d = 2$ align with the numerical results when $d = 1$. That is, the numerical approximation converges with order $s$ for tested initial data in $H^s(2\pi\mathbb{T})$ as long as $1 + \eta^0 > 0$ when using both the smooth and sharp low-pass filters. Figure 5a shows a log-log plot of the relative error measured in the $L^2$-norm, $E_0(U_N)$ and in the $H^1$-norm, $E_1(U_N)$, at time $T = 0.1$ for $2M = 2^j, j = 5, \dots, 9$ and $N = |2M/3|$ with initial data (4.12) with $h_0 = 0.5, u_1 = -v_1 = 0.5, u_h = -v_h = 1$ and $s = 2$. Analogous results also hold for other values of $s$ (we have tested $s = 2.5, 3$).

We have not observed any instabilities in the numerical approximation of system (4.1) due to the use of sharp low-pass filter for any of the initial data we tested satisfying $1 + \eta^0 > 0$. Just as in the case of dimension one, we observe instabilities when violating the non-cavitation assumption $1 + \eta^0 > 0$ for the sharp low-pass filter, but not the smooth low-pass filter, see Figure 6.

For the Hamiltonian Saint-Venant system (4.2), the numerical results are in line with the analysis. Whenever the hyperbolicity condition $1 + \eta^0 - (u^0)^2 - (v^0)^2 > 0$ is satisfied, the numerical scheme converges with the expected rate for both the smooth and the sharp low-pass filter. This is illustrated in Figure 5b for initial data (4.12) with $h_0, u_1 = -v_1 = 0.5, u_h = -v_h = 1$ and $s = 2$.

We observe instabilities in the numerical approximation of system (4.2) when using the sharp low-pass filter for initial data violating $1 + \eta^0 - (u^0)^2 - (v^0)^2 > 0$ but not $1 + \eta^0 > 0$. This is shown in Figure 7. There, we take as initial data (4.12) with $h_0 = 0.5, u_1 = -v_1 = 2, u_h = -v_h = 1$ and $s = 2$. We relate these instabilities to the lack of well-posedness of the underlying system (4.2) when the hyperbolicity condition $1 + \eta^0 - (u^0)^2 - (v^0)^2 > 0$ fails; see Lemma 4.1. The numerical scheme with the smooth low-pass filter does not exhibit instabilities for the tested values of $M$, but limited computational power prevents us from testing very large values of $M$ in dimension two.
<figure><img src="image_11.png" /><figcaption>Figure 5: Plot illustrating the convergence of the numerical schemes (4.8) and (4.7) for the systems (4.1) (in the left) and (4.2) (in the right) in two spatial dimensions as the number of collocation points $2M$ increases. The plot shows the relative error of the numerical solution for initial data (4.12) with $h_0 = 0.5$, $u_1 = -v_1 = 0.5$, $u_h = -v_h = 1$ and $s = 2$ at time $T = 0.1$. The initial data is in $H^2((2\pi\mathbb{T})^3)$ and the relative error is measured in the $L^2$-norm, $E_0$ and in the $H^1$-norm, $E_1$ for $2M = 2^j, j = 6, \dots, 9$ when using either sharp or smooth low-pass filters. To illustrate, the blue and orange lines have slopes $-2$ and $-1$ respectively. The numerical scheme exhibits spectral convergence with both sharp and smooth low-pass filters, for both systems.</figcaption></figure>

## A Technical tools

The following results are standard, and proofs in the Euclidean space (e.g. [14, Theorem 8.3.1]) for product estimates straightforwardly adapt to the periodic setting.

**Proposition A.1 (Continuous embedding).** Let $s \in \mathbb{R}$, $s > d/2$ and $f \in H^s((2\pi\mathbb{T})^d)$. Then $f \in L^\infty((2\pi\mathbb{T})^d)$ and

$$|f|_{L^\infty} \le C(s)|f|_{H^s}.$$

**Proposition A.2 (Interpolation inequality).** Let $s_1, s_2 \in \mathbb{R}$ and $f \in H^{s_1}((2\pi\mathbb{T})^d) \cap H^{s_2}((2\pi\mathbb{T})^d)$. Then for any $0 \le \theta \le 1$, $f \in H^{\theta s_1 + (1-\theta)s_2}((2\pi\mathbb{T})^d)$ and

$$|f|_{H^{\theta s_1 + (1-\theta)s_2}} \le |f|_{H^{s_1}}^{\theta} |f|_{H^{s_2}}^{1-\theta}.$$

**Proposition A.3 (Product estimates).** Let $s_0 > d/2$, $s \ge -s_0$ and $f \in H^s((2\pi\mathbb{T})^d) \cap H^{s_0}((2\pi\mathbb{T})^d)$, $g \in H^s((2\pi\mathbb{T})^d)$. Then $fg \in H^s((2\pi\mathbb{T})^d)$ and

$$|fg|_{H^s} \le C(s_0, s) (|f|_{H^{s_0}} |g|_{H^s} + |f|_{H^s} |g|_{H^{s_0}}).$$

If moreover $s \le s_0$ then

$$|fg|_{H^s} \le C(s_0, s) |f|_{H^{s_0}} |g|_{H^s}.$$

Assuming polynomial nonlinearities, the following proposition is a straightforward consequence of product estimates. Extending this result to general (smooth) functions $P$ requires an analysis that is outside of the scope of the present paper.

**Proposition A.4 (Composition estimates).** Let $s_0 > d/2$, $s \ge -s_0$, $f, g \in H^s((2\pi\mathbb{T})^d) \cap H^{s_0}((2\pi\mathbb{T})^d)$ and $P \in \mathbb{R}[X]$ a polynomial. Then $P(f), P(g) \in H^s((2\pi\mathbb{T})^d)$ and

$$|P(f) - P(g)|_{H^s} \le C(P, s_0, s, |f|_{H^{\max(s_0,s)}}, |g|_{H^{\max(s_0,s)}}) |f - g|_{H^s},$$

$$|P(f) - P(0)|_{H^s} \le C(P, s_0, s, |f|_{H^{s_0}}) |f|_{H^s}.$$
<figure><img src="image_12.png" /><figcaption>Figure 6: Plots of the second derivative $\partial_x^2 u_N$ of the numerical solution to (4.1) for $d=2$ at time $T=0.1$, computed with either the sharp or smooth low-pass filter. The initial data is (4.12) with $s=2$, $u_1 = -v_1 = 0.5$, $u_h = -v_h = 1$, negative minimal depth $h_0 = -0.1$ and $N = \lfloor 2M/3 \rfloor$ for $2M = 2^9$ or $2M = 2^{10}$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 7: Plots of the second derivative $\partial_x^2 u_N$ of the numerical solution to (4.2) for $d = 2$ at time $T = 0.1$, computed with either the sharp or smooth low-pass filter. The initial data is (4.12) with $s = 2$, $u_l = -v_l = 2$, $u_h = -v_h = 1$, positive minimal depth $h_0 = 0.5$ and $N = [2M/3]$ for $2M = 2^9$ or $2M = 2^{10}$.</figcaption></figure>
**Proposition A.5** (Commutator estimates with symbols of order s). Let $s_0 > d/2$, $s \ge 0$ and $\Lambda^s = (\mathrm{Id} - \Delta)^{s/2}$. Let $f \in H^s((2\pi\mathbb{T})^d) \cap H^{s_0+1}((2\pi\mathbb{T})^d)$, $g \in H^{s-1}((2\pi\mathbb{T})^d) \cap H^{s_0}((2\pi\mathbb{T})^d)$. Then

$$|[\Lambda^s, f]g|_{L^2} \le C(s_0, s)(|f|_{H^s} |g|_{H^{s_0}} + |f|_{H^{s_0+1}} |g|_{H^{s_1}}).$$

The following result is shown on $\mathbb{R}^d$ in [8, Lemma 4.5]. The proof straightforwardly adapts to the periodic setting.

**Proposition A.6** (Commutator estimates with operators of order zero). Let $s_0 > d/2$, $s \ge 0$ and $G(D)$ be a Fourier multiplier with symbol $G$ satisfying $|G|_{L^\infty}, ||\cdot||_{\nabla G}|_{L^\infty} \le C_G$. Let $f \in H^{s_0+1}((2\pi\mathbb{T})^d) \cap H^s((2\pi\mathbb{T})^d)$, $g \in H^{s-1}((2\pi\mathbb{T})^d)$. Then

$$|[G(D), f]g|_{H^s} \le C(s_0, s) C_G |f|_{H^{\max(s_0+1,s)}} |g|_{H^{s_1}}.$$

**Remark A.7.** Notice that for smooth symbols considered in this work, namely $S_N(\cdot) = S(\cdot/N)$ where $S$ is even with

$$\begin{cases} S(\mathbf{k}) = 1 & \text{if } \max_{j=1,...,d} |k_j| \le 1/2, \\ S(\mathbf{k}) = 0 & \text{if } \min_{j=1,...,d} |k_j| \ge 1, \\ S(\mathbf{k}) \in [0,1] & \text{otherwise,} \end{cases}$$

and $S^{1/2}$ is Lipschitz-continuous, $S_N^{1/2}$ satisfies the hypotheses of Proposition A.6 uniformly with respect to $N$. Indeed, by Rademacher's theorem we have that $S_N^{1/2}$ is differentiable almost everywhere and its derivative is essentially bounded, and since $S$ has compact support, $|\cdot| \nabla S_N^{1/2} \in L^\infty$. Moreover, we have

$$|S_N^{1/2}|_{L^\infty} + ||\cdot||_{\nabla S_N^{1/2}}|_{L^\infty} = |S^{1/2}|_{L^\infty} + ||\cdot||_{\nabla S^{1/2}}|_{L^\infty}.$$

## Acknowledgements

VD thanks Centre Henri Lebesque ANR-11-LABX-0020-0 for fostering an attractive mathematical environment. JUM acknowledges the support of the project IMod (Grant No. 325114) from the Research Council of Norway.

## References

[1] C. Bardos and E. Tadmor. Stability and spectral convergence of Fourier method for nonlinear problems: on the shortcomings of the 2/3 de-aliasing method. *Numer. Math.*, 129(4):749–782, 2015.

[2] S. Benzoni-Gavage and D. Serre. *Multidimensional hyperbolic partial differential equations. First-order systems and applications*. Oxford Mathematical Monographs. The Clarendon Press Oxford University Press, Oxford, 2007.

[3] M. Cabrera Calvo, F. Rousset, and K. Schratz. Time integrators for dispersive equations in the long wave regime. *Math. Comp.*, 91(337):2197–2214, 2022.

[4] M. Cabrera Calvo and K. Schratz. Uniformly accurate splitting schemes for the Benjamin-Bona-Mahony equation with dispersive parameter. *BIT*, 62(4):1625–1647, 2022.

[5] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. *Spectral methods. Fundamentals in single domains*. Sci. Comput. Berlin: Springer, 2006.

[6] V. A. Dougalis, A. Duran, and L. Saridaki. On the numerical approximation of Boussinesq/Boussinesq systems for internal waves. *Numer. Methods Partial Differential Equations*, 39(5):3677–3704, 2023.

[7] V. Duchêne. Many Models for Water Waves. Open Math Notes, OMN:202109.111309, 2021.
[8] V. Duchêne and B. Melinand. Rectification of a deep water model for surface gravity waves. *Pure Appl. Anal.*, 6(1):73–128, 2024.

[9] V. Duchêne and P. Navaro. Waterwaves1d.jl (version v0.2.1). Zenodo. https://doi.org/10.5281/zenodo.17365679, 2025.

[10] L. Emerald. Local well-posedness result for a class of non-local quasi-linear systems and its application to the justification of Whitham-Boussinesq systems. arXiv preprint:2206.09213.

[11] J. Goodman, T. Hou, and E. Tadmor. On the stability of the unsmoothed Fourier method for hyperbolic equations. *Numer. Math.*, 67(1):93–129, 1994.

[12] D. Gottlieb and J. S. Hesthaven. Spectral methods for hyperbolic problems. volume 128, pages 83–131. 2001. Numerical analysis 2000, Vol. VII, Partial differential equations.

[13] D. Gottlieb and S. A. Orszag. *Numerical analysis of spectral methods: theory and applications*. CBMS-NSF Regional Conference Series in Applied Mathematics, No. 26. Society for Industrial and Applied Mathematics, Philadelphia, PA, 1977.

[14] L. Hörmander. *Lectures on nonlinear hyperbolic differential equations*, volume 26 of *Mathématiques & Applications (Berlin) [Mathematics & Applications]*. Springer-Verlag, Berlin, 1997.

[15] C. Klein, F. Linares, D. Pilod, and J.-C. Saut. On Whitham and related equations. *Stud. Appl. Math.*, 140(2):133–177, 2018.

[16] H.-O. Kreiss and J. Oliger. Comparison of accurate methods for the integration of hyperbolic equations. *Tellus*, 24:199–215, 1972.

[17] D. Lannes. *The water waves problem*, volume 188 of *Mathematical Surveys and Monographs*. American Mathematical Society, Providence, RI, 2013. Mathematical analysis and asymptotics.

[18] A. Majda, J. McDonough, and S. Osher. The Fourier method for nonsmooth initial data. *Math. Comp.*, 32(144):1041–1081, 1978.

[19] G. Métivier. *Para-differential calculus and applications to the Cauchy problem for nonlinear systems*, volume 5 of *Centro di Ricerca Matematica Ennio De Giorgi (CRM) Series*. Edizioni della Normale, Pisa, 2008.

[20] S. A. Orszag. On the elimination of aliasing in finite-difference schemes by filtering high-wavenumber components. *J. Atmos. Sci.*, 28:1074, 1971.

[21] M. O. Paulsen. Long time well-posedness of Whitham-Boussinesq systems. *Nonlinearity*, 35(12):6284–6348, 2022.

[22] J.-C. Saut and L. Xu. The Cauchy problem on large time for surface waves Boussinesq systems. *J. Math. Pures Appl.* (9), 97(6):635–662, 2012.

[23] E. Tadmor. Stability analysis of finite difference, pseudospectral and Fourier-Galerkin approximations for time-dependent problems. *SIAM Rev.*, 29(4):525–555, 1987.

[24] E. Tadmor. Convergence of spectral methods for nonlinear conservation laws. *SIAM J. Numer. Anal.*, 26(1):30–44, 1989.

[25] L. N. Trefethen. *Spectral methods in MATLAB*, volume 10 of *Software, Environments, and Tools*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2000.

[26] J. C. Xavier, M. A. Rincon, D. G. Alfaro Vigo, and D. E. Amundsen. Stability analysis for a fully discrete spectral scheme for Boussinesq systems. *Appl. Anal.*, 97(4):610–632, 2018.