# Stability of finite volume schemes for hyperbolic systems in two space dimensions

Jean-François Coulombel

► To cite this version:

Jean-François Coulombel. Stability of finite volume schemes for hyperbolic systems in two space dimensions. [Research Report] CNRS. 2006. hal-00020015v2

HAL Id: hal-00020015

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

Submitted on 7 Mar 2006

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 CC BY 4.0 - Attribution - International License
# Stability of finite volume schemes for hyperbolic systems in two space dimensions

Jean-François COULOMBEL

CNRS & Université Lille 1, Laboratoire Paul Painlevé, UMR CNRS 8524,

Cité scientifique, 59655 VILLENEUVE D'ASCQ Cedex, France

and Team SIMPAF, INRIA Futurs

Email: jfcoulom@math.univ-lille1.fr

January 11, 2006

## Abstract

We study the stability of finite volume schemes for symmetric hyperbolic systems in two space dimensions, with the Lax-Friedrichs flux. We first show a sufficient condition for the L²-stability of the scheme on a general triangulation. Then we show that this stability condition can be improved when the triangulation is composed of equilateral triangles.

**AMS subject classification:** 65M12, 35L45.

**Keywords:** Hyperbolic systems, finite volume schemes, stability, Lax-Friedrichs flux.

## 1 Introduction, main results

We consider a linear symmetric hyperbolic system in two space dimensions:

$$ egin{cases} 
\partial_t u + A_1 \partial_{x_1} u + A_2 \partial_{x_2} u = 0, & t \ge 0, x \in \mathbb{R}^2, \\ 
u|_{t=0} = u_0, & x \in \mathbb{R}^2, 
\end{cases} 
\quad (1) $$

where $A_1$, and $A_2$ are square $d \times d$ matrices with real coefficients ($d \ge 1$). We assume that $A_1$, and $A_2$ are symmetric, so that the Cauchy problem (1) is well-posed for initial data $u_0 \in L^2(\mathbb{R}^2)$, see e.g. [5]. Moreover, the unique solution to (1) in $C([0, +\infty[; L^2(\mathbb{R}^2))$ satisfies

$$ \forall t \ge 0, \quad \|u(t)\|_{L^2(\mathbb{R}^2)} = \|u_0\|_{L^2(\mathbb{R}^2)}. \quad (2) $$

We introduce a finite volume approximation of (1), with a modified Lax-Friedrichs flux. (We refer to [1, 2, 4] for a general presentation of finite volume schemes for hyperbolic conservation laws.) Let $\mathcal{T}$ denote a triangulation of the plane $\mathbb{R}^2$, that is composed of polygons. For all element $K \in \mathcal{T}$, $\mathcal{N}(K)$ denotes the set of neighbors of $K$. If $L \in \mathcal{N}(K)$, then $\nu_{K,L}$ denotes the outgoing normal vector of $K$ on the edge $K \cap L$. In all what follows, we denote $m(K)$ the area of the polygon $K$, $m(\partial K)$ its perimeter, and $m(K \cap L)$ the length of the edge $K \cap L$ when $L \in \mathcal{N}(K)$. We assume that the triangulation $\mathcal{T}$ satisfies

$$ \sup_{K \in \mathcal{T}} \frac{m(\partial K)}{m(K)} < +\infty. \quad (3) $$

For all vector $\nu = (\nu_1, \nu_2) \in \mathbb{R}^2$, we define the matrix:

$$ A(\nu) := \nu_1 A_1 + \nu_2 A_2. $$
Let $D$ be a positive number, and $\Delta t > 0$ be a time step. Then the finite volume approximation $v$ of the solution $u$ to the Cauchy problem (1) is defined as follows:

$$ \forall n \in \mathbb{N}, \quad \forall K \in \mathcal{T}, \quad \forall (t, x) \in [n\Delta t, (n+1)\Delta t[ \times K, \quad v(t, x) = v_K^n, $$

where, for all $K \in \mathcal{T}$:

$$ \left\{ \begin{aligned} v_K^{n+1} &:= v_K^n - \frac{\Delta t}{m(K)} \sum_{L \in \mathcal{N}(K)} m(K \cap L) \left( \frac{1}{2} A(\nu_{K,L}) (v_K^n + v_L^n) - D(v_L^n - v_K^n) \right), & n \in \mathbb{N}, \\ v_K^0 &:= \frac{1}{m(K)} \int_K u_0(x) dx. \end{aligned} \right. \quad (4) $$

The aim of this work is to derive sufficient conditions on $D$, and $\Delta t$ that ensure the stability of the scheme (4) in the space $L^2(\mathbb{R}^2)$. (In view of the conservation property (2), $L^2(\mathbb{R}^2)$ seems to be the natural framework for this study.) On general triangulations, the way to prove stability is the so-called energy method. This was already used in [6] in the case of the Godunov flux. However, for triangulations that are invariant under sufficiently many translations, the optimal stability criterion follows from Fourier analysis (as for finite difference schemes on cartesian grids). In this case, stability of (4) reduces to the uniform power boundedness of an appropriate symbol (that depends on $D$, $\Delta t$, the triangulation, and on wave numbers $\xi_1, \xi_2$). The derivation, and the analysis of this symbol is the major task of this work.

Let us now state our results. Our first result deals with the case of a general triangulation, that is only assumed to satisfy (3):

**Theorem 1.** Assume that the triangulation $\mathcal{T}$ satisfies (3). Define

$$ \rho_{max} := \sup_{K \in \mathcal{T}} \sup_{L \in \mathcal{N}(K)} \rho(A(\nu_{K,L})), $$

where $\rho$ denotes the spectral radius of a square matrix (with real or complex coefficients). If $\Delta t$, and $D$ satisfy the following inequalities:

$$ \Delta t \rho_{max} \sup_{K \in \mathcal{T}} \frac{m(\partial K)}{m(K)} \le 1, \quad \frac{1}{2} \rho_{max} \le D \le \frac{1}{\Delta t} \inf_{K \in \mathcal{T}} \frac{m(K)}{m(\partial K)} - \frac{1}{2} \rho_{max}, \quad (5) $$

then the finite volume scheme (4) is stable in $L^2(\mathbb{R}^2)$. More precisely, the norm $\|v(t)\|_{L^2(\mathbb{R}^2)}$ is a nonincreasing function of $t$, and

$$ \forall t \ge 0, \quad \|v(t)\|_{L^2(\mathbb{R}^2)} \le \|u_0\|_{L^2(\mathbb{R}^2)}. $$

Observe that the inequalities of (5) make sense because of (3). It is rather remarkable that our condition on the time step $\Delta t$ is the same as the stability condition that was found in [6] for the finite volume scheme with the Godunov flux. The result is surprising since it is known that on cartesian grids, the stability conditions are different for the modified Lax-Friedrichs scheme and for the Godunov scheme. This phenomenon tends to indicate that the energy method does not yield optimal results when the triangulations have symmetries.

Indeed, the limitations (5) can be improved when the triangulation is composed of equilateral triangles, and is oriented as in figure 1. Observe that the triangulation of figure 1 is the only triangulation composed of equilateral triangles (up to rotations and translations). The length of the edges of each equilateral triangle is denoted $h$. For the triangulation of figure 1, we observe that (5) reads:

$$ 4\sqrt{3} \frac{\Delta t}{h} \rho_{max} \le 1, \quad \frac{1}{2} \rho_{max} \le D \le \frac{h}{4\sqrt{3} \Delta t} - \frac{1}{2} \rho_{max}, $$
with

$$ \rho_{max} = \max \left( \rho(A_1), \rho\left(\frac{1}{2}A_1 + \frac{\sqrt{3}}{2}A_2\right), \rho\left(\frac{1}{2}A_1 - \frac{\sqrt{3}}{2}A_2\right) \right). $$

We are going to prove the following result:

**Theorem 2.** Assume that the triangulation $T$ is composed of equilateral triangles, whose edges have length $h$, and is oriented as in figure 1. If $\Delta t$, and $D$ satisfy the following inequalities:

$$ 2\sqrt{3} \frac{\Delta t}{h} \max(\rho(A_1), \rho(A_2)) \le 1, \quad \sqrt{3} \frac{\Delta t}{h} \max(\rho(A_1)^2, \rho(A_2)^2) \le D \le \frac{h}{4\sqrt{3} \Delta t}, \quad (6) $$

then the finite volume scheme (4) is stable in $L^2(\mathbb{R}^2)$. More precisely, the following estimate holds:

$$ \forall t \ge 0, \quad \|v(t)\|_{L^2(\mathbb{R}^2)} \le 2 \|u_0\|_{L^2(\mathbb{R}^2)}. $$

<figure><img src="image_2.png" /><figcaption>Figure 1: Triangulation with equilateral triangles</figcaption></figure>

The rest of this paper is devoted to the proof of Theorems 1, and 2. In section 2, we show Theorem 1 by using the energy method. The analysis closely follows [6]. In section 3, we show Theorem 2. The analysis is divided into several steps. First we compute explicitly the translations that leave the triangulation invariant. We also introduce the decomposition of grid functions into their “even”, and “odd” components. Then we show how the Fourier transform
reduces the investigation of stability to the uniform power boundedness of a suitable symbol. At last, we prove the uniform power boundedness of this symbol under the conditions (6).

In all this article, $\langle x; y \rangle$ denotes the standard hermitian product of two vectors $x, y \in \mathbb{C}^p$, while $\| \cdot \|$ denotes both the associated norm on $\mathbb{C}^p$, and the induced norm on the set of square $p \times p$ matrices (with complex coefficients). The notation is the same for any integer $p$. When vectors have real components, this definition coincides with the standard scalar product, and norm on $\mathbb{R}^p$. Eventually, $I$ denotes the identity matrix (independently of its dimension).

For all integer $n$, $v^n$ is the mesh function that equals $v_K^n$ on the polygon $K$.

## 2 Stability for a general triangulation

We consider the scheme (4) on a triangulation $\mathcal{T}$ that satisfies (3). We also assume that the inequalities (5) are satisfied. To show that the function $(t \mapsto \|v(t)\|_{L^2(\mathbb{R}^2)})$ is nonincreasing, it is sufficient to prove that the sequence $(\|v^n\|_{L^2(\mathbb{R}^2)})_{n \in \mathbb{N}}$ is nonincreasing, or equivalently

$$ \forall n \in \mathbb{N}, \quad \sum_{K \in \mathcal{T}} m(K) \|v_K^{n+1}\|^2 \le \sum_{K \in \mathcal{T}} m(K) \|v_K^n\|^2. $$

Let $n \in \mathbb{N}$, and $K \in \mathcal{T}$. Following the approach of [6], we decompose the vector $v_K^{n+1}$ as the following convex combination:

$$ v_K^{n+1} = \sum_{L \in \mathcal{N}(K)} \frac{m(K \cap L)}{m(\partial K)} v_{K,L}^{n+1}, \quad v_{K,L}^{n+1} := v_K^n + \frac{m(\partial K) \Delta t}{2m(K)} (2DI - A(\nu_{K,L}))(v_L^n - v_K^n). \quad (7) $$

The decomposition (7) can be obtained from (4) by using the well-known relation:

$$ \sum_{L \in \mathcal{N}(K)} m(K \cap L) A(\nu_{K,L}) = A \left( \sum_{L \in \mathcal{N}(K)} m(K \cap L) \nu_{K,L} \right) = 0. \quad (8) $$

Once the vector $v_K^{n+1}$ is written as the convex combination (7), we can use Cauchy-Schwarz' inequality, and derive:

$$ \|v_K^{n+1}\|^2 \le \sum_{L \in \mathcal{N}(K)} \frac{m(K \cap L)}{m(\partial K)} \|v_{K,L}^{n+1}\|^2. \quad (9) $$

It is thus sufficient to estimate the norm of each vector $v_{K,L}^{n+1}$, then (9) will yield an upper bound for $\|v_K^{n+1}\|^2$.

**Estimate of $v_{K,L}^{n+1}$.** Defining the real symmetric matrix

$$ B_{K,L} := \frac{m(\partial K) \Delta t}{2m(K)} [2DI - A(\nu_{K,L})], \quad (10) $$

the vector $v_{K,L}^{n+1}$ reads

$$ v_{K,L}^{n+1} = v_K^n + B_{K,L}(v_L^n - v_K^n). $$

We thus compute:

$$ \begin{aligned} \|v_{K,L}^{n+1}\|^2 - \|v_K^n\|^2 &= \|v_{K,L}^{n+1} - v_K^n\|^2 + 2 \langle v_K^n; v_{K,L}^{n+1} - v_K^n \rangle \\ &= \|B_{K,L}(v_L^n - v_K^n)\|^2 + 2 \langle v_K^n; B_{K,L}(v_L^n - v_K^n) \rangle \\ &= \langle v_L^n - v_K^n; B_{K,L}^2(v_L^n - v_K^n) \rangle - 2 \langle v_K^n; B_{K,L}v_K^n \rangle + 2 \langle v_K^n; B_{K,L}v_L^n \rangle. \end{aligned} $$
We now use the relation

$$2 \langle v_K^n; B_{K,L} v_L^n \rangle = \langle v_K^n; B_{K,L} v_K^n \rangle + \langle v_L^n; B_{K,L} v_L^n \rangle - \langle v_L^n - v_K^n; B_{K,L} (v_L^n - v_K^n) \rangle,$$

and obtain

$$\|v_{K,L}^{n+1}\|^2 - \|v_K^n\|^2 = \langle v_L^n - v_K^n; (B_{K,L}^2 - B_{K,L})(v_L^n - v_K^n) \rangle + \langle v_L^n; B_{K,L}v_L^n \rangle - \langle v_K^n; B_{K,L}v_K^n \rangle. \quad (11)$$

Using the definition (10), and the assumption (5), we get

$$B_{K,L} \geq 0, \quad B_{K,L} - I \leq 0,$$

so the matrix $B_{K,L}^2 - B_{K,L}$ is a nonpositive symmetric matrix. Consequently, (11) yields the following estimate of $v_{K,L}^{n+1}$:

$$\|v_{K,L}^{n+1}\|^2 - \|v_K^n\|^2 \le \langle v_L^n; B_{K,L} v_L^n \rangle - \langle v_K^n; B_{K,L} v_K^n \rangle. \quad (12)$$

**End of the proof of Theorem 1.** We use (12) in the inequality (9), and then multiply by $m(K)$:

$$m(K) \|v_K^{n+1}\|^2 - m(K) \|v_K^n\|^2 \le \sum_{L \in \mathcal{N}(K)} \frac{m(K) m(K \cap L)}{m(\partial K)} \{\langle v_L^n; B_{K,L} v_L^n \rangle - \langle v_K^n; B_{K,L} v_K^n \rangle\}. \quad (13)$$

Using once again the relation (8), and the definition (10) of $B_{K,L}$, we compute

$$\begin{aligned} \sum_{L \in \mathcal{N}(K)} \frac{m(K) m(K \cap L)}{m(\partial K)} \langle v_K^n; B_{K,L} v_K^n \rangle &= \frac{\Delta t}{2} \left\langle v_K^n; \sum_{L \in \mathcal{N}(K)} m(K \cap L) (2D I - A(\nu_{K,L})) v_K^n \right\rangle \\ &= \frac{\Delta t}{2} \langle v_K^n; 2m(\partial K) D v_K^n \rangle \\ &= \frac{\Delta t}{2} \left\langle v_K^n; \sum_{L \in \mathcal{N}(K)} m(K \cap L) (2D I + A(\nu_{K,L})) v_K^n \right\rangle \\ &= \frac{\Delta t}{2} \sum_{L \in \mathcal{N}(K)} m(K \cap L) \langle v_K^n; (2D I - A(\nu_{L,K})) v_K^n \rangle, \end{aligned}$$

where, for the last equality, we have used $\nu_{K,L} = -\nu_{L,K}$. Using this computation in (13), we get

$$\begin{aligned} & m(K) \|v_K^{n+1}\|^2 - m(K) \|v_K^n\|^2 \\ & \le \frac{\Delta t}{2} \sum_{L \in \mathcal{N}(K)} m(K \cap L) \{ \langle v_L^n; (2D I - A(\nu_{K,L})) v_L^n \rangle - \langle v_K^n; (2D I - A(\nu_{L,K})) v_K^n \rangle \}. \quad (14) \end{aligned}$$

Consequently, when we sum the inequalities (14) over the elements $K$ of the triangulation, we end up with

$$\sum_{K \in \mathcal{T}} m(K) \|v_K^{n+1}\|^2 - \sum_{K \in \mathcal{T}} m(K) \|v_K^n\|^2 \le 0.$$

This shows that the function $(t \mapsto \|v(t)\|_{L^2(\mathbb{R}^2)}$ is nonincreasing on $\mathbb{R}^+$. To complete the proof of Theorem 1, we observe that

$$\|v_K^0\|^2 = \frac{1}{m(K)^2} \left\| \int_K u_0(x) dx \right\|^2 \le \frac{1}{m(K)} \int_K \|u_0(x)\|^2 dx,$$

and summing over $K \in \mathcal{T}$, we end up with

$$\|v^0\|_{L^2(\mathbb{R}^2)}^2 = \sum_{K \in \mathcal{T}} m(K) \|v_K^0\|^2 \le \|u_0\|_{L^2(\mathbb{R}^2)}^2. \quad (15)$$

The proof of Theorem 1 is thus complete.
# 3 Stability for a regular triangulation

In this section, we always assume that the triangulation $\mathcal{T}$ is composed of equilateral triangles, whose edges have length $h$, and is oriented as in figure 1. We also assume that the inequalities (6) hold. Observe that the inequality (15) still holds.

## 3.1 A few properties of the triangulation

The first task is to number the triangles as shown in figure 1: for all $(j, k) \in \mathbb{Z}^2$, $\mathcal{T}_{j,k}$ denotes the unique element of $\mathcal{T}$ that is included in the rectangle $[jh\sqrt{3}/2, (j+1)h\sqrt{3}/2] \times [(k-1)h/2, (k+1)h/2]$. The center of the triangle $\mathcal{T}_{j,k}$ is easily computed, and we obtain:

$$ 
\left\{ 
\begin{array}{ll} 
\mathcal{T}_{j,k} = \mathcal{T}_{0,0} + \frac{h}{2}(j\sqrt{3}, k), & \text{if } j+k \text{ is even,} \\ 
\mathcal{T}_{j,k} = \mathcal{T}_{-1,0} + \frac{h}{2}((j+1)\sqrt{3}, k), & \text{if } j+k \text{ is odd.} 
\end{array} 
\right. 
$$

The triangulation $\mathcal{T}$ is thus invariant under the translations

$$ x \in \mathbb{R}^2 \to x + (0, h), \quad \text{and} \quad x \in \mathbb{R}^2 \to x + (h\sqrt{3}, 0). $$

We define:

$$ \mathbb{Z}_e^2 := \{(j,k) \in \mathbb{Z}^2/j+k \text{ is even}\}, \quad \mathbb{Z}_o^2 := \{(j,k) \in \mathbb{Z}^2/j+k \text{ is odd}\}. $$

Let $w$ denote a mesh function, that is a function that is constant on each triangle $\mathcal{T}_{j,k}$, and equals $w_{j,k}$ on $\mathcal{T}_{j,k}$. We adopt the following notation:

$$ \mathbb{E}w(x) := \begin{cases} w_{j,k}, & \text{if } x \in \mathcal{T}_{j,k}, \text{ and } (j,k) \in \mathbb{Z}_e^2, \\ 0, & \text{otherwise,} \end{cases} $$

and

$$ \mathbb{O}w(x) := \begin{cases} w_{j,k}, & \text{if } x \in \mathcal{T}_{j,k}, \text{ and } (j,k) \in \mathbb{Z}_o^2, \\ 0, & \text{otherwise.} \end{cases} $$

When $w \in L^1(\mathbb{R}^2)$, that is, when the sum $\sum_{(j,k)} \|w_{j,k}\|$ is finite, the Fourier transform of $\mathbb{E}w$, and $\mathbb{O}w$ are given by:

$$ \widehat{\mathbb{E}w}(\xi) = \sum_{(j,k) \in \mathbb{Z}_e^2} \exp \left( -ij\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) w_{j,k} \int_{\mathcal{T}_{0,0}} \exp(-i\langle x; \xi \rangle) dx, $$

$$ \widehat{\mathbb{O}w}(\xi) = \sum_{(j,k) \in \mathbb{Z}_o^2} \exp \left( -i(j+1)\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) w_{j,k} \int_{\mathcal{T}_{-1,0}} \exp(-i\langle x; \xi \rangle) dx. $$

Both series are normally convergent for all $\xi \in \mathbb{R}^2$. For later use, we define the following function:

$$ \forall \xi \in \mathbb{R}^2, \quad a(\xi) := \int_{\mathcal{T}_{0,0}} \exp(-i\langle x; \xi \rangle) dx. \qquad (16) $$

Observe that $\mathcal{T}_{-1,0}$ is the symmetric of $\mathcal{T}_{0,0}$ with respect to the origin (see figure 1), so by a change of variables we get

$$ \forall \xi \in \mathbb{R}^2, \quad \int_{\mathcal{T}_{-1,0}} \exp(-i\langle x; \xi \rangle) dx = \overline{a(\xi)}. $$
This simplifies the expression of the Fourier transform of $Ew$, and $Ow$:

$$ 
\begin{aligned}
\widehat{Ew}(\xi) &= a(\xi) \sum_{(j,k) \in \mathbb{Z}_e^2} \exp \left( -ij\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) w_{j,k}, \\
\widehat{Ow}(\xi) &= \overline{a(\xi)} \sum_{(j,k) \in \mathbb{Z}_o^2} \exp \left( -i(j+1)\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) w_{j,k}.
\end{aligned}
\qquad (17) $$

These expressions hold as long as the mesh function $w$ belongs to $L^1(\mathbb{R}^2)$.

When the mesh function $w$ belongs to $L^2(\mathbb{R}^2)$, the functions $Ew$, and $Ow$ are orthogonal in $L^2(\mathbb{R}^2)$, and Plancherel's theorem yields

$$ 
\begin{aligned}
\|w\|_{L^2(\mathbb{R}^2)}^2 &= \|Ew\|_{L^2(\mathbb{R}^2)}^2 + \|Ow\|_{L^2(\mathbb{R}^2)}^2 = \frac{1}{(2\pi)^2} \left( \|\widehat{Ew}\|_{L^2(\mathbb{R}^2)}^2 + \|\widehat{Ow}\|_{L^2(\mathbb{R}^2)}^2 \right) \\
&= \frac{1}{(2\pi)^2} \left\| \left( \begin{matrix} \widehat{Ew} \\ \widehat{Ow} \end{matrix} \right) \right\|_{L^2(\mathbb{R}^2)}^2.
\end{aligned}
$$

Consequently, in order to estimate a mesh function $w$ in $L^2(\mathbb{R}^2)$, we only need to estimate the Fourier transform of its “even”, and “odd” components.

## 3.2 The symbol of the finite volume scheme

For later use, we introduce the notations:

$$ \lambda := \frac{\Delta t}{h}, \quad \alpha := \frac{4}{\sqrt{3}} \lambda D. \qquad (18) $$

From now on, $v_{j,k}^n$ denotes the vector $v_K^n$ when $K$ is the triangle $\mathcal{T}_{j,k}$. With this convention, the finite volume method (4) reads:

$$ 
\left\{ 
\begin{aligned}
v_{j,k}^{n+1} &= v_{j,k}^n - \frac{4\lambda}{\sqrt{3}} \sum_{\mathcal{T}_{l,m} \in \mathcal{N}(\mathcal{T}_{j,k})} \left( \frac{1}{2} A(\nu_{\mathcal{T}_{j,k},\mathcal{T}_{l,m}})(v_{j,k}^n + v_{l,m}^n) - D(v_{l,m}^n - v_{j,k}^n) \right), \quad n \in \mathbb{N}, \\
v_{j,k}^0 &= \frac{4}{h^2\sqrt{3}} \int_{\mathcal{T}_{j,k}} u_0(x) dx.
\end{aligned}
\right. $$

For the triangulation of figure 1, the neighbors of a triangle $\mathcal{T}_{j,k}$, and the outgoing normal vectors are easily computed. For instance, if $j+k$ is even, the neighbors of $\mathcal{T}_{j,k}$ are $\mathcal{T}_{j+1,k}$, $\mathcal{T}_{j,k+1}$, $\mathcal{T}_{j,k-1}$, and the corresponding outgoing normal vectors are $(1,0)$, $(-1,\sqrt{3})/2$, and $(-1,-\sqrt{3})/2$. There are analogous formulae when $(j,k) \in \mathbb{Z}_o^2$. The scheme (4) thus reduces to:

$$ 
\begin{aligned}
v_{j,k}^{n+1} = (1 - 3\alpha)v_{j,k}^n &+ \alpha(v_{j+1,k}^n + v_{j,k-1}^n + v_{j,k+1}^n) \\
& - \frac{2\lambda}{\sqrt{3}} A_1 \left( v_{j+1,k}^n - \frac{1}{2}(v_{j,k-1}^n + v_{j,k+1}^n) \right) - \lambda A_2(v_{j,k+1}^n - v_{j,k-1}^n), \quad \text{if } (j,k) \in \mathbb{Z}_e^2, \quad (19)
\end{aligned}
$$

$$ 
\begin{aligned}
v_{j,k}^{n+1} = (1 - 3\alpha)v_{j,k}^n &+ \alpha(v_{j-1,k}^n + v_{j,k-1}^n + v_{j,k+1}^n) \\
& - \frac{2\lambda}{\sqrt{3}} A_1 \left( \frac{1}{2}(v_{j,k-1}^n + v_{j,k+1}^n) - v_{j-1,k}^n \right) - \lambda A_2(v_{j,k+1}^n - v_{j,k-1}^n), \quad \text{if } (j,k) \in \mathbb{Z}_o^2. \quad (20)
\end{aligned}
$$
Using (19), and (20), we are going to compute the symbol of the finite volume scheme. More precisely, for almost every $\xi \in \mathbb{R}^2$, we want to compute a matrix $\mathcal{G}(\xi)$ that satisfies

$$ \left( \begin{array}{c} \widehat{\mathcal{E}v^{n+1}}(\xi) \\ \widehat{\mathcal{O}v^{n+1}}(\xi) \end{array} \right) = \mathcal{G}(\xi) \left( \begin{array}{c} \widehat{\mathcal{E}v^n}(\xi) \\ \widehat{\mathcal{O}v^n}(\xi) \end{array} \right). $$

Observe that, thanks to (15), we have $v^0 \in L^2(\mathbb{R}^2)$, and a direct induction yields $v^n \in L^2(\mathbb{R}^2)$ for all $n \in \mathbb{N}$. This makes the use of the Fourier transform legitimate.

The existence, and the exact expression of the amplification matrix $\mathcal{G}(\xi)$ will be derived from the following result:

**Lemma 1.** Let $n \in \mathbb{N}$. Then for almost every $\xi = (\xi_1, \xi_2) \in \mathbb{R}^2$, there exists a real number $\vartheta$ so that

$$ \begin{aligned} \widehat{\mathcal{E}v^{n+1}}(\xi) = (1-3\alpha)\widehat{\mathcal{E}v^n}(\xi) &+ \alpha \exp(i\vartheta) \left( \exp(i\xi_1 \frac{h\sqrt{3}}{2}) + 2\cos(\xi_2 \frac{h}{2}) \right) \widehat{\mathcal{O}v^n}(\xi) \\ &- 2\lambda \exp(i\vartheta) \left[ \frac{1}{\sqrt{3}} \left( \exp(i\xi_1 \frac{h\sqrt{3}}{2}) - \cos(\xi_2 \frac{h}{2}) \right) A_1 + i \sin(\xi_2 \frac{h}{2}) A_2 \right] \widehat{\mathcal{O}v^n}(\xi), \quad (21) \end{aligned} $$

and

$$ \begin{aligned} \widehat{\mathcal{O}v^{n+1}}(\xi) = (1-3\alpha)\widehat{\mathcal{O}v^n}(\xi) &+ \alpha \exp(-i\vartheta) \left( \exp(-i\xi_1 \frac{h\sqrt{3}}{2}) + 2\cos(\xi_2 \frac{h}{2}) \right) \widehat{\mathcal{E}v^n}(\xi) \\ &+ 2\lambda \exp(-i\vartheta) \left[ \frac{1}{\sqrt{3}} \left( \exp(-i\xi_1 \frac{h\sqrt{3}}{2}) - \cos(\xi_2 \frac{h}{2}) \right) A_1 - i \sin(\xi_2 \frac{h}{2}) A_2 \right] \widehat{\mathcal{E}v^n}(\xi). \quad (22) \end{aligned} $$

*Proof.* We shall give a detailed proof of (21), and leave the proof of (22) to the interested reader. We assume first of all that the mesh function $v^n$ belongs to $L^1(\mathbb{R}^2) \cap L^2(\mathbb{R}^2)$, and consequently $v^{n+1} \in L^1(\mathbb{R}^2) \cap L^2(\mathbb{R}^2)$. We may thus use (17):

$$ \widehat{\mathcal{E}v^{n+1}}(\xi) = a(\xi) \sum_{(j,k) \in \mathbb{Z}_e^2} \exp \left( -ij\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) v_{j,k}^{n+1}. \quad (23) $$

Using (19) in (23), we first obtain:

$$ \begin{aligned} \widehat{\mathcal{E}v^{n+1}}(\xi) = (1-3\alpha)\widehat{\mathcal{E}v^n}(\xi) &+ \alpha a(\xi) \sum_{(j,k) \in \mathbb{Z}_e^2} \exp \left( -ij\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) (v_{j+1,k}^n + v_{j,k-1}^n + v_{j,k+1}^n) \\ &- \frac{2\lambda}{\sqrt{3}} a(\xi) \sum_{(j,k) \in \mathbb{Z}_e^2} \exp \left( -ij\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) A_1 \left( v_{j+1,k}^n - \frac{1}{2}(v_{j,k-1}^n + v_{j,k+1}^n) \right) \\ &- \lambda a(\xi) \sum_{(j,k) \in \mathbb{Z}_e^2} \exp \left( -ij\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) A_2 (v_{j,k+1}^n - v_{j,k-1}^n). \quad (24) \end{aligned} $$

Observe now that for all $\xi \in \mathbb{R}^2$, there exists a real number $\varphi$ that depends on $\xi$, such that

$$ a(\xi) = \exp(i\varphi) \overline{a(\xi)}. \quad (25) $$
The computation of $\varphi$ seems delicate in the general case, see (16), but the point is that we shall not need the exact expression. Using the relation (25), an intermediate calculation gives the following equalities:

$$
\begin{align}
a(\xi) \sum_{(j,k) \in \mathbb{Z}_e^2} \exp \left( -ij\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) v_{j+1,k}^n & \nonumber \\
&= \exp(i\varphi) \overline{a(\xi)} \sum_{(j,k) \in \mathbb{Z}_e^2} \exp \left( -ij\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) v_{j+1,k}^n & \nonumber \\
&= \exp(i\varphi + i\xi_1 h\sqrt{3}) \overline{a(\xi)} \sum_{(j,k) \in \mathbb{Z}_e^2} \exp \left( -i(j+2)\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) v_{j+1,k}^n & \nonumber \\
&= \exp(i\varphi + i\xi_1 h\sqrt{3}) \overline{a(\xi)} \sum_{(l,k) \in \mathbb{Z}_e^2} \exp \left( -i(l+1)\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) v_{l,k}^n & \nonumber \\
&= \exp(i\varphi + i\xi_1 h\sqrt{3}) \widehat{\mathcal{O}v^n}(\xi). \tag{26}
\end{align}
$$

Similarly, we can obtain the equalities

$$
a(\xi) \sum_{(j,k) \in \mathbb{Z}_e^2} \exp \left( -ij\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) v_{j,k+1}^n = \exp \left( i\varphi + i\xi_1 \frac{h\sqrt{3}}{2} + i\xi_2 \frac{h}{2} \right) \widehat{\mathcal{O}v^n}(\xi), \quad (27)
$$

$$
a(\xi) \sum_{(j,k) \in \mathbb{Z}_e^2} \exp \left( -ij\xi_1 \frac{h\sqrt{3}}{2} - ik\xi_2 \frac{h}{2} \right) v_{j,k-1}^n = \exp \left( i\varphi + i\xi_1 \frac{h\sqrt{3}}{2} - i\xi_2 \frac{h}{2} \right) \widehat{\mathcal{O}v^n}(\xi). \quad (28)
$$

For all $\xi \in \mathbb{R}^2$, we define

$$
\vartheta := \varphi + \xi_1 \frac{h\sqrt{3}}{2}, \tag{28}
$$

with $\varphi$ defined by (25). Then we use (26), and (27) in (24). With the notation (28), we obtain
exactly the relation (21).

When $v^n$ belongs to $L^2(\mathbb{R}^2)$, we approximate $v^n$ by a sequence of mesh functions that belong to $L^1(\mathbb{R}^2) \cap L^2(\mathbb{R}^2)$, and use the continuity of the Fourier transform in $L^2(\mathbb{R}^2)$. $\square$

The equalities (21), (22), can be written in a compact form, up to the introduction of a few
notations:

$$
\begin{align}
z &:= \exp(i\xi_1 \frac{h\sqrt{3}}{2}) + 2\cos(\xi_2 \frac{h}{2}), \nonumber \\
S &:= \frac{1}{\sqrt{3}} \left[ \exp(i\xi_1 \frac{h\sqrt{3}}{2}) - \cos(\xi_2 \frac{h}{2}) \right] A_1 + i \sin(\xi_2 \frac{h}{2}) A_2, \tag{29} \\
\mathcal{G}(\xi) &:= \begin{pmatrix} (1-3\alpha)I & \alpha \exp(i\vartheta)zI \\ \alpha \exp(-i\vartheta)\bar{z}I & (1-3\alpha)I \end{pmatrix} + \begin{pmatrix} 0 & -2\lambda \exp(i\vartheta)S \\ 2\lambda \exp(-i\vartheta)S^* & 0 \end{pmatrix}. \nonumber
\end{align}
$$

The size of the matrix $\mathcal{G}(\xi)$ is $2d \times 2d$. With the notations (29), the relations (21), and (22) can be rewritten as

$$
\left( \begin{array}{c} \widehat{\mathcal{E}v^{n+1}}(\xi) \\ \widehat{\mathcal{O}v^{n+1}}(\xi) \end{array} \right) = \mathcal{G}(\xi) \left( \begin{array}{c} \widehat{\mathcal{E}v^n}(\xi) \\ \widehat{\mathcal{O}v^n}(\xi) \end{array} \right). \qquad (30)
$$

The induction relation (30) immediately implies:

$$
\forall n \in \mathbb{N}, \quad \left( \begin{array}{c} \widehat{\mathcal{E}v^n}(\xi) \\ \widehat{\mathcal{O}v^n}(\xi) \end{array} \right) = \mathcal{G}(\xi)^n \left( \begin{array}{c} \widehat{\mathcal{E}v^0}(\xi) \\ \widehat{\mathcal{O}v^0}(\xi) \end{array} \right).
$$
Observing that Plancherel's Theorem gives

$$2\pi \|v^n\|_{L^2(\mathbb{R}^2)} = \left\| \left( \widehat{\mathbb{E}v^n} \atop \widehat{\mathbb{O}v^n} \right) \right\|_{L^2(\mathbb{R}^2)} = \left\| \mathbb{G}^n \left( \widehat{\mathbb{E}v^0} \atop \widehat{\mathbb{O}v^0} \right) \right\|_{L^2(\mathbb{R}^2)}, \quad (31)$$

we conclude that the $L^2(\mathbb{R}^2)$ stability of the finite volume scheme (4) is equivalent to the uniform power boundedness of the symbol $\mathbb{G}(\xi)$. In the next paragraph, we shall see that the conditions (6) yield the uniform bound 2, which gives

$$2\pi \|v^n\|_{L^2(\mathbb{R}^2)} \le 2 \left\| \left( \widehat{\mathbb{E}v^0} \atop \widehat{\mathbb{O}v^0} \right) \right\|_{L^2(\mathbb{R}^2)} = 2(2\pi) \|v^0\|_{L^2(\mathbb{R}^2)} \le 2(2\pi) \|u_0\|_{L^2(\mathbb{R}^2)},$$

where we have used (15). It is therefore sufficient to derive the uniform bound 2 for $\|\mathbb{G}(\xi)^n\|$ to prove Theorem 2.

### 3.3 Uniform power boundedness of the symbol

The uniform bound for the powers of the symbol will follow from the general result:

**Proposition 1.** Let $T$ be a square $p \times p$ matrix, with complex coefficients, that satisfies

$$\forall X \in \mathbb{C}^p, \quad |\langle TX; X \rangle| \le \|X\|^2. \qquad (32)$$

Then the matrix $T$ satisfies

$$\forall n \in \mathbb{N}, \quad \|T^n\| \le 2.$$

The proof of Proposition 1 may be found in [3]. To end the proof of Theorem 2, we are going to verify that for all $\xi \in \mathbb{R}^2$, the matrix $\mathbb{G}(\xi)$ satisfies the property (32). This verification relies on a simple inequality:

**Lemma 2.** Let $\eta_1, \eta_2, \mu_1, \mu_2 \in \mathbb{R}$, and assume that we have

$$12 \mu_1^2 \le \alpha, \quad 4 \mu_2^2 \le \alpha,$$

where $\alpha > 0$ is defined by (18). Then the following inequality holds:

$$\frac{\alpha}{3} \left\{ (\cos \eta_1 + 2 \cos \eta_2)^2 + \sin^2 \eta_1 \right\} + 8 \mu_1^2 \left\{ (\cos \eta_1 - \cos \eta_2)^2 + \sin^2 \eta_1 \right\} + 8 \mu_2^2 \sin^2 \eta_2 \le 3\alpha.$$

*Proof.* Let us denote

$$I := \frac{\alpha}{3} \left\{ (\cos \eta_1 + 2 \cos \eta_2)^2 + \sin^2 \eta_1 \right\} + 8 \mu_1^2 \left\{ (\cos \eta_1 - \cos \eta_2)^2 + \sin^2 \eta_1 \right\} + 8 \mu_2^2 \sin^2 \eta_2.$$

Then under the assumptions of Lemma 2, we have

$$\begin{align*}
I &= \frac{\alpha}{3} \left\{ 1 + 4 \cos \eta_1 \cos \eta_2 + 4 \cos^2 \eta_2 \right\} + 8 \mu_1^2 \left\{ 1 - 2 \cos \eta_1 \cos \eta_2 + \cos^2 \eta_2 \right\} + 8 \mu_2^2 \sin^2 \eta_2 \\
&= \frac{\alpha}{3} + 8 \mu_1^2 + \left( \frac{4\alpha}{3} - 16 \mu_1^2 \right) \cos \eta_1 \cos \eta_2 + \left( \frac{4\alpha}{3} + 8 \mu_1^2 \right) \cos^2 \eta_1 + 8 \mu_2^2 \sin^2 \eta_2 \\
&\le \frac{\alpha}{3} + 8 \mu_1^2 + \left( \frac{2\alpha}{3} - 8 \mu_1^2 \right) (\cos^2 \eta_1 + \cos^2 \eta_2) + \left( \frac{4\alpha}{3} + 8 \mu_1^2 \right) \cos^2 \eta_1 + 8 \mu_2^2 \sin^2 \eta_2 \\
&= \frac{\alpha}{3} + 8 \mu_1^2 \sin^2 \eta_1 + \frac{2\alpha}{3} \cos^2 \eta_1 + 2\alpha \cos^2 \eta_1 + 8 \mu_2^2 \sin^2 \eta_2 \le \frac{\alpha}{3} + \frac{2\alpha}{3} + 2\alpha = 3\alpha,
\end{align*}$$

which proves the Lemma. □
We note that the inequalities (6) can be rewritten in terms of the numbers $\lambda$, and $\alpha$:

$$4\lambda^2 \rho(A_1)^2 \leq \alpha, \quad 4\lambda^2 \rho(A_2)^2 \leq \alpha, \quad 3\alpha \leq 1. \qquad (33)$$

Let $X \in \mathbb{C}^{2d}$ be a unit vector. We split the vector $X$ in

$$X = \begin{pmatrix} X_e \\ X_o \end{pmatrix}, \quad X_e, X_o \in \mathbb{C}^d, \quad \|X_e\|^2 + \|X_o\|^2 = 1.$$

Using the definition (29), $\mathbb{G}(\xi)$ is decomposed as the sum of a hermitian matrix and a skew-hermitian matrix. We compute

$$\langle \mathbb{G}(\xi)X; X \rangle = 1 - 3\alpha + 2\alpha \operatorname{Re} (z \exp(i\vartheta) \langle X_o; X_e \rangle) - 4i\lambda \operatorname{Im} (\exp(i\vartheta) \langle SX_o; X_e \rangle),$$

which yields

$$|\langle \mathbb{G}(\xi)X; X \rangle|^2 = \left[ 1 - 3\alpha + 3\alpha \frac{2}{3} \operatorname{Re} (z \exp(i\vartheta) \langle X_o; X_e \rangle) \right]^2 + 16\lambda^2 \left[ \operatorname{Im} (\exp(i\vartheta) \langle SX_o; X_e \rangle) \right]^2.$$

Observing that $1 - 3\alpha \geq 0$ because of (33), Cauchy-Schwarz' inequality gives

$$\begin{aligned} |\langle \mathbb{G}(\xi)X; X \rangle|^2 &\leq 1 - 3\alpha + \frac{4\alpha}{3} \left[ \operatorname{Re} (z \exp(i\vartheta) \langle X_o; X_e \rangle) \right]^2 + 16\lambda^2 \|X_e\|^2 \|X_o\|^2 \|S\|^2 \\ &\leq 1 - 3\alpha + \frac{4\alpha}{3} |z|^2 \|X_e\|^2 \|X_o\|^2 + 16\lambda^2 \|X_e\|^2 \|X_o\|^2 \|S\|^2 \\ &\leq 1 - 3\alpha + \frac{\alpha}{3} |z|^2 + 4\lambda^2 \|S\|^2. \end{aligned}$$

Eventually, we use the definition (29) of $z$, and $S$. We have

$$|z|^2 = \left( \cos(\xi_1 \frac{h\sqrt{3}}{2}) + 2 \cos(\xi_2 \frac{h}{2}) \right)^2 + \sin^2(\xi_1 \frac{h\sqrt{3}}{2}),$$

and

$$\begin{aligned} 4\lambda^2 \|S\|^2 &\leq 4 \left( \frac{\lambda}{\sqrt{3}} \left| \exp(i\xi_1 \frac{h\sqrt{3}}{2}) - \cos(\xi_2 \frac{h}{2}) \right| \|A_1\| + \lambda \left| \sin(\xi_2 \frac{h}{2}) \right| \|A_2\| \right)^2 \\ &\leq 8 \mu_1^2 \left| \exp(i\xi_1 \frac{h\sqrt{3}}{2}) - \cos(\xi_2 \frac{h}{2}) \right|^2 + 8 \mu_2^2 \sin^2(\xi_2 \frac{h}{2}), \end{aligned}$$

where we have set

$$\mu_1^2 := \frac{\lambda^2 \|A_1\|^2}{3} = \frac{\lambda^2 \rho(A_1)^2}{3}, \quad \mu_2^2 := \lambda^2 \|A_2\|^2 = \lambda^2 \rho(A_2)^2.$$

Using (33), we have $12\mu_1^2 \leq \alpha$, and $4\mu_2^2 \leq \alpha$. Therefore, we can apply Lemma 2 with $\eta_1 = \xi_1 h \sqrt{3}/2$, and $\eta_2 = \xi_2 h/2$; we obtain

$$\frac{\alpha}{3} |z|^2 + 4\lambda^2 \|S\|^2 \leq 3\alpha.$$

This shows that the symbol $\mathbb{G}(\xi)$ satisfies the property (32). Thanks to Proposition 1, the norm $\|\mathbb{G}(\xi)^n\|$ is bounded by 2 for all $n$, and for all $\xi$. Using this bound in (31), and then applying Plancherel's Theorem, we can complete the proof of Theorem 2.
**Acknowledgments** Research of the author was supported by the EU financed network HYKE, HPRN-CT-2002-00282. The author warmly thanks Catalin Badea for pointing out the result of Proposition 1, and for stimulating discussions.

**References**

[1] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In *Handbook of numerical analysis, Vol. VII*, pages 713–1020. North-Holland.

[2] E. Godlewski, P.-A. Raviart. *Numerical approximation of hyperbolic systems of conservation laws*. Springer-Verlag, 1996.

[3] M. Goldberg, E. Tadmor. On the numerical radius and its applications. *Linear Algebra Appl.*, 42:263–284, 1982.

[4] D. Kröner. *Numerical schemes for conservation laws*. John Wiley & Sons Ltd., 1997.

[5] A. Majda. *Compressible fluid flow and systems of conservation laws in several space variables*. Springer-Verlag, 1984.

[6] J.-P. Vila, P. Villedieu. Convergence of an explicit finite volume scheme for first order symmetric systems. *Numer. Math.*, 94(3):573–602, 2003.