Rge Nic Se

C Aise

Ina P Rist

Ignot T I
email: pignotti@univaq.it

Università Di L'aquila

Dipartimento Di

Ingegneria E Scienze

A Note on the Asymptotic Stability of Wave-Type Equations with Switching Time-Delay

Keywords: some concrete models are described Mathematics Subject Classification 2000. 35L05, 93D15 Wave equation, stabilization, switching time delay type feedbacks

We consider second-order evolution equations with intermittently delayed/ not-delayed damping. We give sufficient conditions for asymptotic and exponential stability, completing our previous results from Nicaise and

Introduction

Let H be a real Hilbert space and let A : D(A) → H be a positive self-adjoint operator with a compact inverse in H. Denote by V := D(A 1 2 ) the domain of A 1 2 . Let B 1 , B 2 be time-dependent linear operators, B 1 ∈ L(H, H) and B 2 (t) : U → H, where U is a real Hilbert space with norm and inner product denoted respectively by ∥ • ∥ U and ⟨•, •⟩ U . We assume that B 1 and B 2 act alternately, that is

B * 1 (t)B * 2 (t) = 0, ∀t > 0.
Let us consider the problem u tt (t) + Au(t) + B 1 (t)B * 1 (t)u t (t) + B 2 (t)B * 2 (t)u t (tτ ) = 0 t > 0 (6.1) u(0) = u 0 and u t (0) = u 1 (6.2) where the constant τ > 0 is the time delay. Time delay effects appear in many applications and practical problems and it is by now well known that even an arbitrarily small delay in the feedback may destabilize a system which is uniformly exponentially stable in absence of delay. For some examples of this destabilizing effect of time delays we refer to [START_REF] Datko | Not all feedback stabilized hyperbolic systems are robust with respect to small time delays in their feedbacks[END_REF][START_REF] Datko | An example on the effect of time delays in boundary feedback stabilization of wave equations[END_REF][START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF][START_REF] Xu | Stabilization of wave systems with input delay in the boundary control[END_REF].

In [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF] we considered the wave equation with both dampings acting simultaneously, that is B 1 (t) = µ 1 and B 2 (t) = µ 2 , with µ 1 , µ 2 ∈ IR + , and we proved that if µ 1 > µ 2 then the system is uniformly exponentially stable. Otherwise, if µ 2 ≥ µ 1 , that is the delay term prevails on the not-delayed one, then there are instability phenomena, namely, there are unstable solutions for arbitrarily small (large) delays.

The stabilization problem for second-order evolution equations with switching time delay is, in some sense, related to the one for second-order evolution equations damped by positive/negative feedbacks (see [START_REF] Haraux | Asymptotic stability for intermittently controlled second-order evolution equations[END_REF]). See also [START_REF] Pignotti | A note on stabilization of locally damped wave equations with time delay[END_REF] for the relation between wave equation with time delay in the damping and wave equation with indefinite damping, i.e., damping which changes sign in different subsets of the domain.

We firstly studied this subject in [START_REF] Nicaise | Asymptotic stability of second-order evolution equations with intermittent delay[END_REF], in a more general setting with respect to the problem considered here. Assuming that an observability inequality holds for the conservative model associated with (6.1) and (6.2) and, through the definition of a suitable energy, we obtained sufficient conditions ensuring asymptotic stability. Under more restrictive assumptions exponential stability estimates were also obtained.

An analogous problem has been considered in [START_REF] Ammari | Stabilization by switching time-delay[END_REF] for 1D models for the wave equation but with a different approach. Indeed, in [START_REF] Ammari | Stabilization by switching time-delay[END_REF] we obtain stability results for particular values of the time delays, related to the length of the domain, by using the D'Alembert formula.

The results of [START_REF] Nicaise | Asymptotic stability of second-order evolution equations with intermittent delay[END_REF] are improved in [START_REF] Nicaise | Stability results for second-order evolution equations with switching time-delay[END_REF] by removing a quite restrictive assumption on the size of the "bad" terms, i.e., the terms with time delay (cf. assumption (3.3) of [START_REF] Nicaise | Asymptotic stability of second-order evolution equations with intermittent delay[END_REF]). Indeed, expected from [START_REF] Haraux | Asymptotic stability for intermittently controlled second-order evolution equations[END_REF] and from the relation between delay problems and problems with antidamping, the delay feedback operator B 2 may be also large but then has to act on small time intervals. Moreover, in [START_REF] Nicaise | Stability results for second-order evolution equations with switching time-delay[END_REF] we consider also the case when B 1 is unbounded. In this short note, we complete the results from [START_REF] Nicaise | Asymptotic stability of second-order evolution equations with intermittent delay[END_REF][START_REF] Nicaise | Stability results for second-order evolution equations with switching time-delay[END_REF] by improving the sufficient conditions that guarantee stability results but by staying in the case of bounded operators B 1 .

This chapter is organized as follows. In Section 6.2, we recall the wellposedness result proved in [START_REF] Nicaise | Asymptotic stability of second-order evolution equations with intermittent delay[END_REF] and the observability estimate in short-time proved by Haraux, Martinez, and Vancostenoble. In Section 6.3, we obtain asymptotic and exponential stability results for the abstract model under suitable conditions, completing the results of [START_REF] Nicaise | Stability results for second-order evolution equations with switching time-delay[END_REF]. Finally, in Section 6.4, we give some concrete applications of our results.

Well-Posedness

In this section we recall a well-posedness result, proved by the authors in [START_REF] Nicaise | Asymptotic stability of second-order evolution equations with intermittent delay[END_REF], for problem (6.1) and (6.2).

We assume that for all n ∈ IN, there exists t n > 0 with t n < t n+1 and such that

B 2 (t) = 0 ∀t ∈ I 2n = [t 2n , t 2n+1 ), B 1 (t) = 0 ∀t ∈ I 2n+1 = [t 2n+1 , t 2n+2 ), with B 1 ∈ C 1 ([t 2n , t 2n+1 ]; L(H, H)) and B 2 ∈ C([t 2n+1 , t 2n+2 ]; L(U, H)).
Moreover, we assume that τ ≤ T 2n for all n ∈ IN, where T n denotes the length of the interval I n , i.e.,

T n = t n+1 -t n , n ∈ IN . (6.3)
Under these assumptions, the following result holds (see Theorem 2.1 in [START_REF] Nicaise | Asymptotic stability of second-order evolution equations with intermittent delay[END_REF]). Theorem 6.1 Under the above assumptions, for any u 0 ∈ V and u 1 ∈ H, the system (6.1) and (6.2)

has a unique solution u ∈ C([0, ∞); V) ∩ C 1 ([0, ∞); H).
To obtain stability results for problem (6.1) and (6.2), we assume now that for all n ∈ IN, there exist three positive constants m 2n , M 2n , and M 2n+1 with m 2n ≤ M 2n and such that for all u ∈ H we have

(i) m 2n ∥u∥ 2 H ≤ ∥B * 1 (t)u∥ 2 H ≤ M 2n ∥u∥ 2 H for t ∈ I 2n = [t 2n , t 2n+1 ), ∀n ∈ IN; (ii) ∥B * 2 (t)u∥ 2 U ≤ M 2n+1 ∥u∥ 2 H for t ∈ I 2n+1 = [t 2n+1 , t 2n+2 ), ∀n ∈ IN.
Denote by E S (•) the standard energy for wave-type equations, i.e.,

E S (t) = E S (u; t) := 1 2 ( ∥u(t)∥ 2 V + ∥u t (t)∥ 2 H
) .

Observe that, on the time intervals I 2n , n ∈ IN, only the standard dissipative damping acts. So, the following observability estimate holds (see theorem 3.1 of [START_REF] Haraux | Asymptotic stability for intermittently controlled second-order evolution equations[END_REF]).

Proposition 6.2 Assume (i).

There exists a constant c, independent of the length T 2n of the interval I 2n , such that for any solution of (6.1) and ( 6.2) it holds

E S (t 2n+1 ) ≤ 1 1 + c m 2n T -3 2n +T -1 2n +M 2n m 2n T -1 2n E S (t 2n ), n ∈ IN. (6.4)

Stability Results

Now, as in [START_REF] Nicaise | Stability results for second-order evolution equations with switching time-delay[END_REF], we assume

inf n∈IN m 2n M 2n+1 > 0. (6.5) Note that assumption (3.3) in [10] is instead equivalent to inf n∈IN m 2n M 2n+1 > 1 2 .
Let us introduce the energy of the system

E(t) = E(u; t) := 1 2 ( ∥u(t)∥ 2 V + ∥u t (t)∥ 2 H + ξ 2 ˆt t-τ ∥B * 2 (s + τ )u t (s)∥ 2 U ds
) , (6.6) where ξ is a positive number satisfying

ξ < inf n∈IN m 2n M 2n+1 . (6.7)
The following estimates are proved as in proposition 3.1 of [START_REF] Nicaise | Stability results for second-order evolution equations with switching time-delay[END_REF].

Proposition 6.3 Assume (i), (ii), and (6.5). For any regular solution of problem (6.1) and (6.2) the energy is decreasing on the intervals I 2n , n ∈ IN, and

E ′ (t) ≤ - m 2n 2 ∥u t ∥ 2 H . (6.8)
Moreover, on the intervals

I 2n+1 , n ∈ IN, E ′ (t) ≤ M 2n+1 2 ( ξ + 1 ξ ) ∥u t ∥ 2 H . (6.9)
The following theorem, due to the more general assumption (6.5) and due to the more general setting, improves and generalizes theorem 3.3 of [START_REF] Nicaise | Asymptotic stability of second-order evolution equations with intermittent delay[END_REF]. Theorem 6.4 Assume (i), (ii), and (6.5).

If ∞ ∑ n=0 M 2n+1 T 2n+1 < +∞, (6.10
)

and ∞ ∑ n=0 ln   1 1 + c m 2n T -3 2n +T -1 2n +M 2n m 2n T -1 2n + τ ξM 2n+1   = -∞, (6.11)
then system (6.1) and (6.2) is asymptotically stable, that is any solution u of (6.1) and (6.2) satisfies E S (u, t) → 0 for t → +∞.

Proof Note that (6.9) implies

E ′ (t) ≤ M 2n+1 ( ξ + 1 ξ ) E(t), t ∈ I 2n+1 = [t 2n+1 , t 2n+2 ), n ∈ IN.
Then we have

E(t 2n+2 ) ≤ e (ξ+ 1 ξ )M 2n+1 T 2n+1 E(t 2n+1 ), ∀n ∈ IN. (6.12)
Now, note that

E(t 2n+1 ) = E S (t 2n+1 ) + ξ 2 ˆt2n+1 t 2n+1 -τ ∥B * 2 (s + τ )u t (s)∥ 2 U ds.
Now observe that, since T 2n ≥ τ, n ∈ IN, the variable s + τ in the above integral belongs to

I 2n+1 ∪ I 2n+2 . If s + τ ∈ I 2n+1 , then ∥B * 2 (s + τ )u t (s)∥ 2 U ≤ M 2n+1 ∥u t (s)∥ 2 H . Otherwise, if s + τ ∈ I 2n+2 , then B * 2 (s + τ )u t (s) = 0. Then, E(t 2n+1 ) ≤ E S (t 2n+1 ) + ξ 2 M 2n+1 ˆt2n+1 t 2n+1 -τ ∥u t (s)∥ 2 H ds. (6.13)
Now, since (t 2n+1τ, t 2n+1 ) ⊂ I 2n and in I 2n the system is dissipative, from (6.13) we deduce

E(t 2n+1 ) ≤ E S (t 2n+1 ) + τ M 2n+1 ξE S (t 2n+1 -τ ) ≤ E S (t 2n+1 ) + τ M 2n+1 ξE S (t 2n ). (6.14)
Combining Proposition 6.2 and (6.14) we then obtain

E(t 2n+1 ) ≤   1 1 + c m 2n T -3 2n +T -1 2n +M 2n m 2n T -1 2n + τ ξM 2n+1   E S (t 2n ), (6.15)
and therefore

E S (t 2n+2 ) ≤ E(t 2n+2 ) ≤ e (ξ+ 1 ξ )M 2n+1 T 2n+1   1 1 + c m 2n T -3 2n +T -1 2n +M 2n m 2n T -1 2n + τ ξM 2n+1   E (t 2n ). (6.16)
Since (6.16) holds for any n ∈ IN we can deduce

E S (t 2n+2 ) ≤ n ∏ p=0 e (ξ+ 1 ξ )M2p+1T2p+1   1 1 + c m2p T -3 2p +T -1 2p +M2pm2pT -1 2p + τ ξM 2p+1   E S (0). ( 6
E S (t) ≤ E(t) ≤ e (ξ+ 1 ξ )M 2n+1 T 2n+1 E(t 2n+1 ), (6.19)
where in the second inequality we have used (6.12). Then, by (6.15) and (6.17)-(6.19), asymptotic stability occurs if

∞ ∑ n=0 [( ξ + 1 ξ ) M 2n+1 T 2n+1 + ln cn ] = -∞, (6.20)
where

cn = 1 1 + c m 2n T -3 2n +T -1 2n +M 2n m 2n T -1 2n + τ ξM 2n+1 . (6.21) Thus, if ∞ ∑ n=0 M 2n+1 T 2n+1 < +∞, ∞ ∑ n=0 ln cn = -∞,
system (6.1) and (6.2) is asymptotically stable.

We now show that under additional assumptions on the coefficients T n , m n , M n an exponential stability result holds. Theorem 6.5 Assume (i), (ii), and (6.5). Assume also that where cn is as in (6.21). Then, there exist two positive constants γ, µ such that E S (t) ≤ γe -µt E S (0), t > 0, (6.25)

T 2n = T * ∀n ∈ IN, ( 6
for every solution of problem (6.1) -(6.2).

Proof From (6.16) and ( 6.24) we obtain

E S (T * + T) ≤ dE S (0),
and also

E S (n(T * + T)) ≤ d n E S (0), ∀n ∈ IN.
Then, the standard energy E S (•) satisfies an exponential estimate like (6.25) (see lemma 1 of [START_REF] Gugat | Boundary feedback stabilization by time delay for one-dimensional wave equations[END_REF]).

Remark 6.6

In the assumptions of Theorem 6.5, from (6.17) we can see that exponential stability also holds if instead of (6.24) we assume

∃n ∈ IN such that k(n+1)+n ∏ p=k(n+1)
e (ξ+ 1 ξ )M 2p+1 Tc p ≤ d < 1, ∀k = 0, 1, 2, . . .

Examples

Here, we illustrate some concrete examples falling in our previous abstract setting, namely the wave equation, the elasticity system, the Midlin-Timoshenko model, the Petrowsky system.

The Wave Equation

As

u tt (x, t) -∆u(x, t) + b 1 (t)u t (x, t) + b 2 (t)χ ω u t (x, t -τ ) = 0 in Ω × (0, +∞) (6.26) u(x, t) = 0 on ∂Ω × (0, +∞) (6.27) u(x, 0) = u 0 (x) and u t (x, 0) = u 1 (x) in Ω (6.28) with initial data (u 0 , u 1 ) ∈ H 1 0 (Ω) × L 2 (Ω) and b 1 , b 2 in L ∞ (0, +∞) such that b 1 (t)b 2 (t) = 0, ∀t > 0.
Moreover, we assume

(i w ) 0 < m 2n ≤ b 1 (t) ≤ M 2n , b 2 (t) = 0, for all t ∈ I 2n = [t 2n , t 2n+1 ), and b 1 ∈ C 1 ( Ī2n ), for all n ∈ IN; (ii w ) |b 2 (t)| ≤ M 2n+1 , b 1 (t) = 0, for all t ∈ I 2n+1 = [t 2n+1 , t 2n+2 ), and b 2 ∈ C( Ī2n+1 ), for all n ∈ IN.
This problem enters into our previous framework, if we take H = L 2 (Ω) and the operator A defined by

A : D(A) → H : u → -∆u, where D(A) = H 1 0 (Ω) ∩ H 2 (Ω).
The operator A is a self-adjoint and positive operator with a compact inverse in H and is such that V = D(A 1/2 ) = H 1 0 (Ω). We then define the operator B 1 as

B 1 : H → H : v → √ b 1 v, (6.29)
and, denote U := L 2 (Ω), the operator B 2 as

B 2 : U → H : v → √ b 2 ṽ, (6.30)
where ṽ ∈ L 2 (Ω) is the extension of v by zero outside ω.

It is easy to verify that B 1 B * 1 (φ) = b 1 φ and B 2 B * 2 (φ) = b 2 φχ ω , for φ ∈ H. This shows that problem (6.26) and (6.28) enters in the abstract framework (6.1) and (6.2). Moreover, (i w ) and (ii w ) easily imply (i) and (ii) of Section 6.3. Therefore we can restate Proposition 6.3. Now, the energy functional is

E(t) = 1 2 ˆΩ{u 2 t (x, t) + |∇u(x, t)| 2 }dx + ξ 2 ˆt t-τ |b 2 (s + τ )| ˆωu 2 t (x,
E ′ (t) ≤ - m 2n 2 ˆΩ u 2 t (x, t)dx. (6.32)
Moreover, on the intervals

I 2n+1 , n ∈ IN, E ′ (t) ≤ M 2n+1 2 ( ξ + 1 ξ )ˆΩ u 2 t (x, t)dx. ( 6.33)
Thus, the stability results of Theorems 6.4 and 6.5 apply to this model.

The Elasticity System

In the same setting than in the previous section, we consider the following elastodynamic system b 1 , b 2 satisfying the same assumptions as in Section 6.4.1. Here the state variable u is vector-valued and λ, µ are the Lamé coefficients that are positive real numbers. Finally for a (smooth enough) vector-valued function

u tt (x, t) -µ∆u(x, t) -(λ + µ)∇ div u(x, t) + b 1 (t)u t (x, t) + b 2 (t)χ ω u t (x, t -τ ) = 0 in Ω × (0, +∞) (6.34) u(x, t) = 0 on ∂Ω × (0, +∞) (6.35) u(x, 0) = u 0 (x) and u t (x, 0) = u 1 (x) in Ω (6.36) with initial data (u 0 , u 1 ) ∈ H 1 0 (Ω) d × L 2 (Ω) d and
v : Ω → R d , div v is its standard divergence, namely div v = d ∑ j=1 ∂ j v j .
As before this problem enters into our abstract setting, once we take H = L 2 (Ω) d , and A defined by

A : D(A) → H : u → -µ∆u(x, t) -(λ + µ)∇ div u, where D(A) = H 1 0 (Ω) d ∩ H 2 (Ω) d .
The operator A is a self-adjoint and positive operator with a compact inverse in H and is such that

V = D(A 1/2 ) = H 1 0 (Ω) d equipped with the inner product (u, v) V = ˆΩ ( µ d ∑ i,j=1 ∂ i u j ∂ i v j + (λ + µ) div u div v ) dx, ∀u, v ∈ H 1 0 (Ω) d .
We further define U = L 2 (ω) d and the operators B i , i = 1, 2, as follows:

B 1 : U → H : v → √ b 1 v, B 2 : U → H : v → √ b 2 ṽ,
ṽ being again the extension of v by zero outside ω. As before

B 1 B * 1 (φ) = b 1 φ, and B 2 B * 2 (φ) = b 2 φχ ω ,
for any φ ∈ H. So, problem (6.34)-( 6.36) enters in the abstract framework (6.1) and (6.2). Moreover, (i w ) and (ii w ) easily imply (i) and (ii) of Section 6.2. Therefore, the results of Section 6.3 apply also to the system (6.34)-(6.36) with the energy defined by

E(t) = 1 2 ˆΩ { |u t | 2 (x, t) + µ d ∑ i,j=1 (∂ i u j (x, t)) 2 + (λ + µ)(div u(x, t)) 2 } dx + ξ 2 ˆt t-τ |b 2 (s + τ )| ˆω |u t (x, s)| 2 dx ds.

The Mindlin-Timoshenko Model

In the same setting as in Section 6.4.1, we consider the internal stabilization of the following Mindlin-Timoshenko (beam/plate) model (for similar models, see chapter 5 of [START_REF] Lagnese | Boundary stabilization of thin plates[END_REF], chapters 2 and 4 of [START_REF] Lagnese | Modelling analysis and control of thin plates[END_REF], [START_REF] Fernández Sare | On the stability of Mindlin-Timoshenko plates[END_REF][START_REF] Nicaise | Internal stabilization of a Mindlin-Timoshenko model by interior feedbacks[END_REF]).

w tt (x, t) = div(K(∇w + u))(x, t) -a 1 (t)w t (x, t)
-a 2 (t)χ ω w t (x, tτ ) in Ω × (0, +∞) (6.37)

u tt (x, t) = div Cϵ(u)(x, t) -K(∇w + u)(x, t) -b 1 (t)u t (x, t) -b 2 (t)χ ω u t (x, t -τ ) in Ω × (0, +∞) (6.38)
u(x, t) = 0, w(x, t) = 0 on ∂Ω × (0, +∞) (6.39) u(x, 0) = u 0 (x), u t (x, 0) = u 1 (x), w(x, 0) = w 0 (x), w t (x, 0) = w 1 (x) in Ω.

(6.40)

If d = 1 (resp. d = 2
) the scalar variable w represents the displacement of the beam (resp. plate) in the vertical direction, while the vectorial variable u = (u i ) d i=1 is the angles of rotation of a filament of the beam (resp. plate). Here K belongs to L ∞ (Ω) d×d , is symmetric and positive definite, i.e., there exists a positive constant k 0 such that

X ⊤ K(x)X ≥ k 0 , ∀X ∈ IR d , for a.e. x ∈ Ω, while C = (c ijkℓ ) is a tensor such that c ijkℓ = c jiℓk = c kℓij ∈ L ∞ (Ω), (6.41)
all indices running over the integers 1, . . . , d. These quantities are related to the constitutive materials of the beam/plate. As usual for u = (u i ) d i=1 , ϵ(u) is the linear strain tensor defined by As usual we assume that C is positive definite in the sense that there exists µ 0 > 0 such that C(x)ϵ : ϵ ≥ µ 0 |ϵ| 2 , ∀ϵ ∈ IR d×d , for a.e. x ∈ Ω.

ϵ(u) = (ϵ ij (u)) d i,j=1 with ϵ ij (u) = 1 2 (∂ i u j + ∂ j u i ).
(6.42)

We further recall that for a (smooth enough) matrix-valued function w = (w ij ) : Ω → R d×d , div w is its divergence line by line, i.e., div w =

( d ∑ j=1 ∂ j w ij ) d i=1 .
Finally we require that the functions b 1 , b 2 satisfy the assumptions (i w ) and (ii w ) from Section 6.4.1, and similarly for a 1 and a 2 , we suppose that (i mt ) 0 < m 2n ≤ a 1 (t) ≤ M 2n , a 2 (t) = 0, for all t ∈ I 2n = [t 2n , t 2n+1 ), and a 1 ∈ C 1 ( Ī2n ), for all n ∈ IN; ṽ being again the extension of v by zero outside ω. Clearly one has B 1 B * 1 (u, w) = (b 1 u, a 1 w) and B 2 B * 2 (u, w) = (b 2 uχ ω , a 2 wχ ω ), for any (u, w) ∈ H. So, problem (6.37)-(6.39) enters in the abstract framework (6.1) and (6.2). Moreover, (i w ), (ii w ), (i mt ), and (ii mt ) imply (i) and (ii) from Section 6.2, and consequently the results of Section 6.3 apply also to this system with the energy defined by The operator A is self-adjoint and positive, has a compact inverse in H and satisfies D(A 1/2 ) = H 2 (Ω) ∩ H 1 0 (Ω). We then define U = L 2 (ω) and the operators B i , i = 1, 2, by (6.29) and (6.30). So, problem (6.43)-(6.45) enters in the abstract framework (6.1) and (6.2). Moreover, (i w ) and (ii w ) easily imply (i) and (ii) of Section 6.2.

E(t) = 1
Therefore, the results of Section 6.3 apply also to the plate model.

  .22) with T * ≥ τ, and T 2n+1 = T ∀n ∈ IN. (6.23)Moreover, assume thatsup n∈IN e (ξ+ 1 ξ )M 2n+1 Tc n = d < 1,(6.24)

For a d

 d × d matrix ϵ = (ϵ ij ) d i,j=1 the product Cϵ = ((Cϵ) ij ) d i,j=1 is the d × d matrix given by (Cϵ) ij = d ∑ k,ℓ=1 c ijkℓ ϵ kℓ .

  u t (x, tτ ) = 0 in Ω × (0, +∞) (6.43) u(x, t) = ∆u(x, t) = 0 on ∂Ω × (0, +∞) (6.44) u(x, 0) = u 0 (x) and u t (x, 0) = u 1 (x) in Ω (6.45) with initial data (u 0 , u 1 ) ∈ H 2 (Ω) ∩ H 1 0 (Ω) × L 2 (Ω)and b 1 , b 2 satisfying the same assumptions as in Section 6.4.1. Now, we take H = L 2 (Ω) and let A be the operatorA : D(A) → H : u → ∆ 2 u, (6.46) where D(A) = {v ∈ H 1 0 (Ω) ∩ H 4 (Ω) : ∆u = 0 on ∂Ω}.

  a concrete application let us consider the wave equation with internal damping. More precisely, let Ω ⊂ IR d (with d a positive natural number) be an open bounded domain with a boundary ∂Ω of class C 2 . We denote by ω a subset of Ω.Let us consider the initial boundary value problem

  Assume (i w ), (ii w ), and (6.5). Then, for every regular solution of problem (6.26)-(6.28) the energy is decreasing on the intervals I 2n , n ∈ IN, and

	s)dx ds.
	(6.31)
	Proposition 6.7

4 The Petrovsky System Let

  (s + τ )||u t (x, s)| 2 + |a 2 (s + τ )||w t (x, s)| 2 )dx ds. Ω ⊂ IR d bean open bounded set with a boundary ∂Ω of class C 4 and let ω be any fixed subset of Ω .Let us consider the initial boundary value problem u

	2	ˆΩ{ (Cε(u)(x, t)) : ε(ū(x, t)))+K(∇w+u)(x, t)) • (∇w+ū)(x, t)) }	dx
	+	ξ 2 ˆt t-τ	ˆω(|b 2

6.4.tt (x, t) + ∆ 2 u(x, t) + b 1 (t)u t (x, t) + b 2 (t)χ ω

(ii mt ) |a 2 (t)| ≤ M 2n+1 , a 1 (t) = 0, for all t ∈ I 2n+1 = [t 2n+1 , t 2n+2 ), and a 2 ∈ C( Ī2n+1 ), for all n ∈ IN.

System (6.37)-(6.39) can be viewed as a coupling between the wave equation in w with the dynamical elastic system in u with internal feedbacks with delays.

Again this problem enters into our abstract setting by using Friedrichs extension theorem. Namely we take H = L 2 (Ω) d × L 2 (Ω) with its natural inner product

that is clearly compactly embedded into H and the sesquilinear and symmetric form

with U = (u, w), U * = (u * , w * ) ∈ V. Indeed using Korn's and Poincaré's inequalities, it is not difficult to check that this sesquilinear form is coercive in V, namely there exists α > 0 such that

Hence it is well known that the operator A associated with the triple (a, V, H) is a self-adjoint and positive operator with a compact inverse in H with D(A 1/2 ) = V. This operator is defined by

and

By the definition of a, it is easy to see that

and then

We further define U = H and the operators B i , i = 1, 2, as follows:


Christophe Berthon
email: christophe.berthon@univ-nantes.fr

Solène Bulteau
email: solene.bulteau1@gmail.com

Françoise Foucher
email: francoise.foucher@ec-nantes.fr

Meissa M'baye
email: meissaths@gmail.com

Michel-Dansac Victor
email: victor.michel-dansac@inria.fr

¶‖

A very easy high-order well-balanced reconstruction for hyperbolic systems with source terms

When adopting high-order finite volume schemes based on MUSCL reconstruction techniques to approximate the weak solutions of hyperbolic systems with source terms, the preservation of the steady states turns out to be very challenging. Indeed, the designed reconstruction must preserve the steady states under consideration in order to get the required well-balancedness property. A priori, to capture such a steady state, one needs to solve some strongly nonlinear equations. Here, we design a very easy correction to high-order finite volume methods. This correction can be applied to any scheme of order greater than or equal to 2, such as a MUSCL-type scheme, and ensures that this scheme exactly preserves the steady solutions. The main discrepancy with usual techniques lies in avoiding the inversion of the nonlinear function that governs the steady solutions. Moreover, for under-determined steady solutions, several nonlinear functions must be considered simultaneously. Since the derived correction only considers the evaluation of the governing nonlinear functions, we are able to deal with under-determined stationary systems. Several numerical experiments illustrate the relevance of the proposed well-balanced correction, as well as its main limitation, namely the fact that it may fail at being both well-balanced and more than second-order accurate for a specific class of initial conditions.

Introduction 1.General framework

The present work is devoted to the numerical approximation of the weak solutions of an evolution law of the form ∂ t w + ∂ x f (w) = S(w, x), x ∈ R, t > 0, (

where w : R × R + → Ω ⊂ R N denotes the unknown vector. The set Ω stands for the set of the admissible states and it is assumed to be convex. The flux function f : Ω → R n is assumed to be smooth enough, say C 1 . The source term S : Ω × R → R N is also supposed to be smooth enough for w ∈ Int Ω. We also emphasize the need for the function x ∈ R → S(w, x) to be a known smooth function for all w ∈ Ω. For stability reasons, in the present work, the matrix ∇ w f (w) is assumed to be diagonalizable in R so that the homogeneous system extracted from (1.1) is hyperbolic. The PDE system (1.1) is endowed with initial data w(x, t = 0) = w 0 (x), where w 0 (x) ∈ Ω for all x ∈ R is a given function.

Because of the source term S(w, x), there exists non-trivial steady solutions governed by

∂ x f (w) = S(w, x), w(x 0 ) = w 0 , (1.2)
where w 0 is a given state in Ω and x 0 a given point in R. Now, if the above system can be integrated, there exists G : Ω × R → R N such that the stationary solutions are governed by ∂ x G(w, x) = 0, w(x 0 ) = w 0 .

In fact, it is not always possible to integrate (1.2). Even then, G is not necessarily unique and, generally, it must be restricted according to some invariant domains. Usually, the steady solutions are restricted to some particular definition of G. Then, we have to deal with a sequence (G ℓ ) 1≤ℓ≤L with L < +∞ given.

Equipped with these comments, in this work, we only consider steady solutions defined by

∂ x G ℓ (w, x) = 0, 1 ≤ ℓ ≤ L, w(x 0 ) = Π eq ℓ (w 0 ), (1.3)
where we have denoted by Π eq (w) the projection of w over the invariant domain under consideration. More broadly, each steady state is characterized by two sets of equations. The first one, involving G, controls the space variations of the steady-state solution. The second one, which involves Π, controls the solution in a pointwise fashion.

Illustrating models

In order to illustrate the relevance of such a definition of the steady states, let us present some examples of particular interest. First, let us adopt the well-known shallow water model given by

w = h q , f (w) =   q q 2 h + g h 2 2   , S(w, x) = 0 -gh∂ x z , (1.4)
where g > 0 is the gravity constant, z(x) the given smooth topography function, h is the water height and q is the water discharge. The smooth steady solutions (see [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF]) are given by

     ∂ x q = 0, ∂ x q 2 2h 2 + g(h + z) = 0, (1.5)
with w(x 0 ) = w 0 for a given w 0 in Ω, where the set Ω of admissible states is defined as follows:

Ω = t (h, q) ∈ R 2 ; h ≥ 0, q ∈ R . (1.6)
As a consequence, we immediately obtain G(w, x) =   q q 2 2h 2 + g(h + z)   and Π eq (w) = w.

(1.7)

Now, if we restrict the definition of the steady states by only considering the usual lake at rest given by q = 0,

∂ x (h + z) = 0,
then we easily get G(w, x) = 0 h + z and Π eq (w) = h 0 .

(1.8)

The second example we present is given by the shallow water equations with a Manning friction source term (see [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography or Manning friction[END_REF][START_REF] Manning | On the flow of water in open channels and pipes[END_REF][START_REF] Delestre | SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies[END_REF]) and a flat bottom. This model reads

w = h q , f (w) =   q q 2 h + g h 2 2   , S(w, x) =   0 -κ q|q| h η   , (1.9)
where κ > 0 is the friction coefficient and η ̸ = 1 is the Manning exponent. The set of the admissible states is given by (1.6). After [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography or Manning friction[END_REF], the stationary solutions are given as follows:

     ∂ x q = 0, ∂ x -q 2 h η-1 η -1 + g h η+2 η + 2 + κxq|q| = 0,
with w(x 0 ) = w 0 for a given w 0 ∈ Ω, so that we immediately obtain

G(w, x) =    q -q 2 h η-1 η -1 + g h η+2 η + 2 + κxq|q|  
 and Π eq (w) = w.

(1.10)

Next, we present a system involving a non-unique definition of G, the Euler model with gravity (see [START_REF] Desveaux | A well-balanced scheme to capture non-explicit steady states in the Euler equations with gravity[END_REF][START_REF] Xu | A well-balanced kinetic scheme for gas dynamic equations under gravitational field[END_REF]). This model reads as follows:

w =    ρ q E    , f (w) =       q q 2 ρ + p (E + p) q ρ       , S(w, x) =    0 -ρ∂ x φ -q∂ x φ    , (1.11)
where φ : R → R stands for a given smooth gravitational potential and p := p(ρ, E -1 2 q 2 ρ ) denotes the pressure law, with E the total energy. The set of admissible states is defined here by Ω = t (ρ, q, E) ∈ R 3 ; ρ > 0, q ∈ R, E -1 2

q 2 ρ > 0 .
Concerning the steady solutions, we are concerned by steady solutions at rest governed by (see [START_REF] Desveaux | A well-balanced scheme to capture non-explicit steady states in the Euler equations with gravity[END_REF][START_REF] Xu | A well-balanced kinetic scheme for gas dynamic equations under gravitational field[END_REF])

q = 0, ∂ x p + ρ∂ x φ = 0, with w(x 0 ) =    ρ 0 0 E 0    .
Once again, the system to govern the steady-state solutions turns out to be under-determined and we have to focus on particular families of steady solutions. According to [START_REF] Desveaux | A well-balanced scheme to capture non-explicit steady states in the Euler equations with gravity[END_REF][START_REF] Xu | A well-balanced kinetic scheme for gas dynamic equations under gravitational field[END_REF], three families of steady states are of prime importance. The first family is given by

     q = 0, ∂ x ρ = 0, ∂ x (p + ρφ) = 0, with w(x 0 ) =    ρ 0 0 E 0    ,
for all ρ 0 > 0 and E 0 > 0.

In order to define both second and third families of steady states, we have to impose that the function E → p(ρ, E) is invertible and we denote by p -1 E (ρ, .) this inverse function so that p -1 E (ρ, p(ρ, E)) = E. Now, the second steady state family reads

     q = 0, ∂ x (p -κρ) = 0, ∂ x (φ + κlnρ) = 0, with w(x 0 ) =    ρ 0 0 p -1 E (ρ 0 , κρ 0 )    ,
where κ > 0 is a given parameter. The last steady state family is defined by

         q = 0, ∂ x (p -κρ γ ) = 0, ∂ x κγ γ -1 ρ γ-1 + φ = 0, with w(x 0 ) =   ρ 0 0 p -1 E (ρ 0 , κρ γ 0 )   , (1.12)
where γ > 1 is a given parameter.

As a consequence, we get

G 1 (w, x) =    q ρ p + ρφ    , G 2 (w, x) =    q p -κρ φ + κlnρ    , G 3 (w, x) =      q p -κρ γ κγ γ -1 ρ γ-1 + φ      , (1.13) with Π eq 1 (w) =    ρ 0 E    , Π eq 2 (w) =    ρ 0 p -1 E (ρ, κρ)    and Π eq 3 (w) =    ρ 0 p -1 E (ρ, κρ γ )    . (1.14)

Main motivation

Now, considering the derivation of numerical schemes approximating the solutions of (1.1) and able to accurately, or even exactly, capture the steady solutions defined by (1.3) has been an important challenge during the two last decades. Numerous techniques have been designed for the shallow water model with topography (1.4) supplemented by the lake at rest (1.8). For a non-exhaustive bibliography, the reader is referred to [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF][START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Parés | On the well-balance property of Roe's method for nonconservative hyperbolic systems. Applications to shallow-water systems[END_REF][START_REF] Gallardo | On a well-balanced high-order finite volume scheme for shallow water equations with topography and dry areas[END_REF][START_REF] Berthon | A positive preserving high order VFRoe scheme for shallow water equations: a class of relaxation schemes[END_REF][START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF][START_REF] Bryson | Well-balanced positivity preserving centralupwind scheme on triangular grids for the Saint-Venant system[END_REF][START_REF] Berthon | Efficient well-balanced hydrostatic upwind schemes for shallow-water equations[END_REF][START_REF] Chen | A new hydrostatic reconstruction scheme based on subcell reconstructions[END_REF][START_REF] Cheng | A new approach for designing moving-water equilibria preserving schemes for the shallow water equations[END_REF]. More recently, in [START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF][START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[END_REF][START_REF] Xing | On the advantage of well-balanced schemes for moving-water equilibria of the shallow water equations[END_REF][START_REF] Xing | Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium[END_REF][START_REF] Berthon | A fully well-balanced, positive and entropy-satisfying godunov-type method for the shallow-water equations[END_REF][START_REF] Berthon | Fully well-balanced, positive and simple approximate Riemann solver for shallow water equations[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography or Manning friction[END_REF], extensions are given in order to deal with moving steady states given by (1.7). In [START_REF] Bouchut | Upwinding of the source term at interfaces for Euler equations with high friction[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography or Manning friction[END_REF][START_REF] Guermond | Well-balanced second-order finite element approximation of the shallow water equations with friction[END_REF], the Manning-type friction source term is adopted and suitable discretizations are introduced to capture steady states given by (1.10). Regarding the discretization of the Euler model with gravity (1.11), the reader is referred to [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF][START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter-dependent source. The case of Euler system with friction[END_REF][START_REF] Xu | A well-balanced kinetic scheme for gas dynamic equations under gravitational field[END_REF][START_REF] Luo | A well-balanced symplecticity-preserving gas-kinetic scheme for hydrodynamic equations under gravitational field[END_REF][START_REF] Xing | High order well-balanced WENO scheme for the gas dynamics equations under gravitational fields[END_REF][START_REF] Käppeli | Well-balanced schemes for the Euler equations with gravitation[END_REF][START_REF] Desveaux | A well-balanced scheme to capture non-explicit steady states in the Euler equations with gravity[END_REF][START_REF] Berberich | Second Order Finite Volume Scheme for Euler Equations with Gravity which is Well-Balanced for General Equations of State and Grid Systems[END_REF][START_REF] Li | Well-balanced discontinuous Galerkin methods with hydrostatic reconstruction for the Euler equations with gravitation[END_REF][START_REF] Thomann | A second-order positivity-preserving well-balanced finite volume scheme for Euler equations with gravity for arbitrary hydrostatic equilibria[END_REF] where numerical strategies are developed to capture steady states according to the pairs (G ℓ , Π eq ℓ ) defined by (1.13) -(1.14). In the present work, we are not interested in the derivation of well-balanced schemes, namely schemes able to capture steady solutions given by (1.3). Here, our purpose concerns the high-order extensions obtained by involving a polynomial reconstruction procedure. Indeed, as soon as the well-balancedness property must be preserved, the reconstruction may involve strong difficulties. In particular, to be well-balanced, the usual reconstruction approaches need to invert G ℓ (w, x) with respect to w for one given ℓ as long as the function w → G ℓ (w, x) is invertible. We immediately understand that dealing simultaneously with distinct functions G ℓ (w, x) does not seem reachable. Moreover, imposing that the application w → G ℓ (w, x) is invertible is a strong assumption, not satisfied in general by physical models.

In this work, we present a very easy strategy to force any reconstruction procedure to preserve the steady solutions given by (1.3) just evaluating the applications G ℓ (w, x) according to the projection Π eq ℓ (w). Let us emphasize that this technique still requires the prior knowledge of a well-balanced first-order scheme, which may be a challenging endeavor. To address such an issue, the present work is organized as follows. In order to set the framework and the main notations, in Section 2 we introduce the numerical schemes and the usual MUSCL second-order strategy [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF]. In addition, we present the main issues when enforcing the polynomial reconstruction to be well-balanced. Section 3 is then devoted to the strategy designed here, which ensures that the expected well-balancedness property is satisfied by any reconstruction. The proposed improvement comes from a suitable evaluation of the pairs (G ℓ , Π eq ℓ ) 1≤ℓ≤L ; G ℓ is never inverted. In Section 4, we present a high-order well-balanced extension. Finally, Section 5 is devoted to several numerical experiments to illustrate the relevance of the designed high-order reconstruction improvement. We also illustrate the current main limitation of the scheme. Indeed, starting from a perturbed steady state solution, we have so far not been able to correctly discretize the initial condition to achieve both the well-balance property and an order of accuracy greater than two. Section 6 concludes this discussion and suggests some perspectives.

Issues of the well-balanced second-order MUSCL schemes

To approximate the solutions of (1.1), the space is discretized by introducing a sequence of cells (x i-1 2 , x i+ 1 2 ), for all i ∈ Z, with a constant size ∆x. We denote by

x i = (x i-1 2 + x i+ 1
2 )/2 the center of each cell. We set t n+1 = t n + ∆t to discretize the time domain with a time step ∆t. In general, ∆t is restricted according to a CFL condition.

At time t n , we denote by w n i a constant approximation of the solution of (1.1) over the cell (x i-1 2 , x i+ 1 2 ). To evolve this approximation in time, we adopt a finite volume scheme of the form

w n+1 i = w n i - ∆t ∆x f n i+ 1 2 -f n i-1 2 + ∆t 2 S n i-1 2 + S n i+ 1 2 , (2.1)
where we have set

f n i+ 1 2 = f ∆ (w n i , w n i+1 ) and S n i+ 1 2 = S ∆ (w n i , w n i+1 , x i , x i+1 , ∆x).
In order to get a consistent scheme, the numerical flux function f ∆ and the discrete source term S ∆ are assumed to be Lipschitz-continuous and to verify f ∆ (w, w) = f (w) and S ∆ (w, w, x, x, 0) = S(w, x).

(2.2)

At this level, the reader is referred to the large literature devoted to the derivation of a well-balanced scheme according to the system of interest. Here, we have imposed the well-balancedness property according to the definition (1.3) of the steady states. As a consequence, we get w n+1 i = w n i as long as, for all i in Z, we have

G ℓ (w n i , x i ) = G ℓ (w n i+1 , x i+1 ) and w n i = Π eq ℓ (w n i ) with 1 ≤ ℓ ≤ L, (2.3)
to enforce the invariant domain according to (1.3). Now, we focus on a second-order extension, see for instance [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF][START_REF] Berthon | An entropy preserving MOOD scheme for the Euler equations[END_REF]. To address such an issue, we have to introduce suitable reconstructed states, denoted by w ± i+ 1 2 , on each side of the interface x i+ 1 2 . This reconstruction is said to be second-order accurate in space if, for all i ∈ Z, we have

w - i+ 1 2 = w(x i+ 1 2 , t n ) + O(∆x 2 ) and w + i+ 1 2 = w(x i+ 1 2 , t n ) + O(∆x 2 ), (2.4)
for some smooth function x → w(x, t n ) such that

w n i = 1 ∆x x i+ 1 2 x i-1 2 w(x, t n )dx.
Equipped with this second-order reconstruction, from the first-order scheme (2.1), we define a second-order scheme as follows:

w n+1 i = w n i - ∆t ∆x f ± i+ 1 2 -f ± i-1 2 + ∆tS ± i , (2.5)
where we have set f ±

i+ 1 2 = f ∆ (w - i+ 1 2 , w + i+ 1 2
), and where S ± i is a second-order approximation of the source term average, i.e.

S ± i = 1 ∆x x i+ 1 2 x i-1 2 S(w(x, t n ), x) dx + O(∆x 2 ). (2.6)
A classical choice for such second-order accurate schemes is to use the second-order midpoint approximation:

S ± i = S 1 2 w + i-1 2 + w - i+ 1 2 , x i . (2.7)
It is clear that second-order accuracy is achieved as soon as the reconstructed states are defined. At the interface x i+ 1 2 , the reconstructed states read (for instance, see [START_REF] Berthon | An entropy preserving MOOD scheme for the Euler equations[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF])

w - i+ 1 2 = w n i + 1 2 L(w n i -w n i-1 , w n i+1 -w n i ), w + i+ 1 2 = w n i+1 - 1 2 L(w n i+1 -w n i , w n i+2 -w n i+1 ), (2.8)
where L : R N × R N → R N are Lipschitz-continuous functions, which satisfy

L(w, w) = w for all w ∈ R N , ∃M > 0 such that ∥L(w L , w R )∥ ≤ M max(∥w L ∥, ∥w R ∥), ∀ w L , w R ∈ R N .
A large body of literature is devoted to introduce suitable definition of L (for instance, see [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] and references therein). Now, by adopting (2.8), the steady states are, in general, not preserved by such a reconstruction. Indeed, in order to recover the expected well-balancedness property, we require

w - i+ 1 2 = w n i and w + i+ 1 2 = w n i+1 for all i ∈ Z, (2.9)
as soon as (w n i ) i∈Z defines a steady state according to (2.3). Except for linear steady states, the steady condition (2.9) is lost. As a consequence, a particular attention must be paid on the definition of L to preserves the steady states.

Currently, the reconstruction on G ℓ instead of w is preferred (see [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[END_REF]). For a fixed ℓ, denoted by ℓ ⋆ , we perform the reconstruction as follows:

G - i+ 1 2 = G ℓ ⋆ (w n i , x i ) + 1 2 L G ℓ ⋆ (w n i , x i ) -G ℓ ⋆ (w n i-1 , x i-1 ), G ℓ ⋆ (w n i+1 , x i+1 ) -G ℓ ⋆ (w n i , x i ) , G + i+ 1 2 = G ℓ ⋆ (w n i+1 , x i+1 ) - 1 2 L G ℓ ⋆ (w n i+1 , x i+1 ) -G ℓ ⋆ (w n i , x i ), G ℓ ⋆ (w n i+2 , x i+2 ) -G ℓ ⋆ (w n i+1 , x i+1 ) ,
The reconstructed states w ± i+ 1 2 at the interface x i+ 1 2 are then defined by

   G ℓ ⋆ (w - i+ 1 2 , x i+ 1 2 ) = G - i+ 1 2 , G ℓ ⋆ (w + i+ 1 2 , x i+ 1 2 ) = G + i+ 1 2 .
(2.10)

We immediately remark that (2.9) holds as soon as (w n i ) i∈Z defines a steady state for G ℓ ⋆ according to (2.3). However, the function w → G ℓ ⋆ (w, x) must be inverted. Such a procedure may be very costly, or even impossible to carry out if G ℓ ⋆ (., x) is not invertible.

In fact, in the simpler situation of the lake at rest for the shallow water equation, where G is given by (1.8), we have to solve a linear 2 × 2 system. But for a moving steady state, i.e. with G defined by (1.7), the uniqueness of the reconstructed states w ± i+ 1 2 is no longer ensured. Next, considering (1.10), neither the existence nor the uniqueness of w ±

i+ 1 2 is ensured.
Moreover, adopting such a procedure needs to fix ℓ. As a consequence, it is not possible to deal with steady states governed by several families G ℓ (w, x) with 1 ≤ ℓ ≤ L for L ≥ 2. Such a restriction arises for instance for the Euler equations with gravity, where we consider three steady state families.

To summarize the failure of the usual well-balanced reconstruction technique, the reconstructed states, solution of (2.10), may not exist or not be unique. Moreover, since we have to solve a nonlinear system, the evaluation of the reconstructed states turns out to be computationally expensive. In addition, such a reconstruction technique preserves only one steady state family while some systems involve several steady state families.

A very easy well-balanced MUSCL reconstruction

The objective is now to derive a reconstruction technique able to preserve the steady states but never involving an inversion of G ℓ . To address such an issue, we suggest to improve the usual reconstruction (2.8) as follows:

w-

i+ 1 2 = w n i + 1 2 θ n i+ 1 2 L(w n i -w n i-1 , w n i+1 -w n i ),
w+

i+ 1 2 = w n i+1 - 1 2 θ n i+ 1 2 L(w n i+1 -w n i , w n i+2 -w n i+1 ), (3.1)
where the correction θ n i+ 1 2 must be an approximation of 1, with at least second-order accuracy, which vanishes for pairs (w n i , w n i+1 ) satisfying (2.3). We propose the following formulation of θ n i+ 1 2 :

θ n i+ 1 2 = ε n i+ 1 2 ε n i+ 1 2 + ∆x C n i+ 1 2 k , with (3.2)
ε n i+ 1 2 = L ℓ=1 ∥G ℓ (w n i+1 , x i+1 ) -G ℓ (w n i , x i )∥ + ∥w n i+1 -Π eq ℓ (w n i+1 )∥ + ∥w n i -Π eq ℓ (w n i )∥ , (3.3)
where k ≥ 2 must be selected and where C n i+ 1 2 ̸ = 0 is any expression independent from ∆x. We shall suggest an expression of C n i+ 1 2 in the numerical experiments. From now on, it is worth noting that ε n i+ 1 2 = 0 if and only if the pair (w n i , w n i+1 ) defines a local steady state, according to (2.3), at the interface x i+ 1 2 . Concerning the source term discretization, we adopt the following definition:

S± i = 1 2 1 -θ n i-1 2 S n i-1 2 + 1 -θ n i+ 1 2 S n i+ 1 2 + 1 2 θ n i-1 2 + θ n i+ 1 2 S ± i , (3.4)
where S ± i is given by (2.7) and where S n i± 1 2 comes from the first-order discretization (2.1). As a consequence, the second-order MUSCL scheme now reads

w n+1 i = w n i - ∆t ∆x f ± i+ 1 2 -f ± i-1 2 + ∆t S± i , (3.5)
where we have set

f ± i+ 1 2 = f ∆ ( w- i+ 1 2 , w+ i+ 1 2
).

(3.6)

Before we establish the main properties satisfied by the second-order MUSCL scheme (3.5) -(3.6) with the reconstructed states (3.1) and the source term discretization (3.4), let us recall the definition of the order of accuracy that is adopted here (for instance, see [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]). Definition 1. For some smooth solution w(x, t) of (1.1), let us consider

w n i = 1 ∆x x i+ 1 2 x i-1 2 w(x, t n ) dx. (3.7)
Define w n+1 i by (2.5). The scheme (2.5) is said to be of order τ in time and δ in space if, for all i in Z, we have

w n+1 i = 1 ∆x x i+ 1 2 x i-1 2 w(x, t n + ∆t) dx - ∆t ∆x F i+ 1 2 -F i-1 2 + ∆tS i , (3.8)
where

F i+ 1 2 = O(∆t τ ) + O(∆x δ ) and S i = O(∆t τ ) + O(∆x δ
). Now, arguing the above definition of the order of accuracy, the improved reconstruction technique based on θ i+ 1 2 is established to yield a second-order accurate and well-balanced scheme. Theorem 2. Given a well-balanced first-order scheme, the scheme (3.5) -(3.6), with reconstructed states given by (3.1) and a source term discretization given by (3.4), satisfies the following properties:

(i) it is second-order accurate in space for unsteady solutions;

(ii) it is well-balanced, i.e. it exactly preserves steady solutions: if (w n i ) i∈Z defines a steady state according to (2.3), then w n+1 i = w n i for all i in Z;

(iii) it is robust, i.e. if the original reconstruction (2.8) preserves Ω, then Ω remains invariant by the improved reconstruction (3.1).

Proof. We establish properties (i), (ii) and (iii) in order.

(i) To establish the order of accuracy, let us consider w(x, t) a smooth unsteady solution of (1.1). By integration of (1.1) over

(x i-1 2 , x i+ 1 2 ) × (t n , t n + ∆t), we get 1 ∆x x i+ 1 2 x i-1 2 w(x, t n + ∆t) dx - 1 ∆x x i+ 1 2 x i-1 2 w(x, t n ) dx + ∆t ∆x 1 ∆t t n +∆t t n f w(x i+ 1 2 , t) dt - 1 ∆t t n +∆t t n f w(x i-1 2 , t) dt = ∆t 1 ∆t∆x t n +∆t t n x i+ 1 2 x i-1 2
S(w(x, t), x) dx dt. does not vanish, we have θ n

F i+ 1 2 = f ± i+ 1 2 - 1 ∆t ∆t 0 f w(x i+ 1 2 , t n + t) dt, S i = S± i - 1 ∆t∆x ∆t 0 x i+ 1 2 x i-1 2 S(w(x, t n + t), x) dx dt, where f ±
i+ 1 2 = 1 + O(∆x k ).
As a consequence, in the current unsteady context, we get w-

i+ 1 2 = w - i+ 1 2
+ O(∆x k ) and w+

i+ 1 2 = w + i+ 1 2 + O(∆x k ),
where w ± i+ 1 2 are given by (2.8), and with k ≥ 2. Since (2.4) holds for the second-order polynomial reconstruction, we immediately obtain w-

i+ 1 2 = w(x i+ 1 2 , t n ) + O(∆x 2
) and w+

i+ 1 2 = w(x i+ 1 2 , t n ) + O(∆x 2 ).
Assuming a Lipschitz-continuous numerical flux function such that the consistency condition (2.2) holds, we have = 0 for all i in Z. As a consequence, we have θ n i+ 1 2 = 0, which leads to w-

f ± i+ 1 2 = f (w(x i+ 1 2 , t n )) + O(∆x
i+ 1 2 = w n i and w+ i+ 1 2 = w n i+1 , while S± i = 1 2 (S n i-1 2 + S n i+ 1 2
). Put in other words, the reconstruction vanishes for steady states. Then, the original well-balanced first-order scheme (2.1) is recovered and the steady states are preserved.

(iii) We finally turn to the robustness of the improved reconstructed states w±

i+ 1 2
. We remark that w-

i+ 1 2 = (1 -θ n i+ 1 2 )w n i + θ n i+ 1 2 w - i+ 1 2
and w+

i+ 1 2 = (1 -θ n i+ 1 2 )w n i+1 + θ n i+ 1 2 w + i+ 1 2 , where θ n i+ 1 2
, defined by (3.2), belongs to [0, 1], and where w ±

i+ 1 2
are given by the initial reconstruction (2.8).

Since the states w n i , w n i+1 and w ± i+ 1 2 belong to Ω, the states w±

i+ 1 2
turn out to be convex combinations of states in Ω. With Ω a convex set, we immediately deduce that w±

i+ 1 2 are in Ω.
The proof is thus completed.

To conclude this section, we emphasize that we have designed a well-balanced reconstruction procedure by only evaluating (G ℓ (w n i , x i )) 1≤ℓ≤L and never solving some nonlinear equations. Moreover, the introduced procedure simultaneously deals with all the involved families of steady states and it is not necessary to give more importance to one than to another.

Well-balanced high-order finite volume extension

The above well-balanced improvement for the second-order MUSCL scheme is easily extended to yield a well-balanced and high-order accurate scheme. Note that in the high-order case, particular attention should be paid to the initialization. Indeed, on the one hand, using usual high-order initialization with cell averages will render the scheme unable to exactly preserve steady states given pointwise by (1.3). On the other hand, using pointwise initialization will cap the scheme to second-order accuracy. Instead, we propose the following way to compute the discrete initial condition in cell i:

w 0 i = (1 -θ i )w 0 (x i ) + θ i 1 ∆x x i+1/2 x i-1/2 w 0 (x) dx, (4.1)
where

θ i = max(θ i-1/2 , θ i+1/2
). This initial condition is given by a cell average far from a steady state and by pointwise values close to a steady state, which lifts both the issues mentioned above.

Remark 3. We remark that, in some situations, even this initialization procedure may fail to be high-order accurate and accurate. For instance, consider an initial condition made of a perturbation of a steady solution.

In some cell i where the perturbation starts, one interface (for instance x i-1 2 ) corresponds to a steady solution, while the other one (for instance x i+ 1 2 ) corresponds to the unsteady perturbation. If θ i = 0, then the high-order approximation of the perturbation at interface x i+ 1 2 is lost; otherwise, the steady solution at interface x i-1 2 is lost. Also, if θ i = 0, the second-order accurate zone will extend by one cell every time iteration. This means that, due to this failure in initialization that we have, so far, been unable to correct, the well-balanced scheme may be limited to second-order accuracy in experiments consisting in perturbed steady solutions. This situation is reflected in the upcoming Theorem 4 and in the numerical experiments from Section 5.2.4.

Now, to build the well-balanced and high-order accurate scheme, let us first introduce a high-order reconstruction according to existing work, see for instance [START_REF] Diot | Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF][START_REF] Diot | The multidimensional optimal order detection method in the threedimensional case: very high-order finite volume method for hyperbolic systems[END_REF]. With w(x, t) a given smooth function, we define w n i by adopting (3.7). Now, we consider the following polynomial reconstruction of degree d in space:

p n w (x; i) = w n i + π w i (x -x i ), (4.2)
where π w i is a polynomial function of degree d such that, for all x ∈ (x i-1 2 , x i+ 1 2 ), we have

p n w (x; i) = w(x, t n ) + O(∆x d+1
) and 1 ∆x

x i+ 1 2 x i-1 2 p n w (x; i) dx = w n i . (4.3)
Equipped with this reconstruction of degree d, a scheme of space order δ = d + 1 is derived. The reader is referred to [START_REF] Diot | Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF][START_REF] Diot | The multidimensional optimal order detection method in the threedimensional case: very high-order finite volume method for hyperbolic systems[END_REF] where such reconstruction techniques are derived.

From this high-order reconstruction, we now give the associated high-order well-balanced scheme to approximate the weak solutions of (1.1) as follows:

w n+1 i = w n i - ∆t ∆x f ± i+ 1 2 -f ± i-1 2 + ∆t S± i , (4.4)
with the numerical flux function given by (3.6), with

w - i+ 1 2 = (p n w ) + i and w + i+ 1 2 = (p n w ) - i+1 , (4.5)
where (p n w ) ± i is the following well-balanced modification of the high-order polynomial reconstruction (4.2) evaluated at the interface point x i± 1 2 :

(p n w ) ± i = w n i + θ n i± 1 2 π w i ± ∆x 2 , (4.6)
with θ n i+ 1 2 defined by (3.2) -(3.3). In (3.2), the parameter k must be fixed larger than δ = d + 1 in order to preserve the order δ of the polynomial reconstruction. Concerning the source term approximation, we start with an approximation of order δ of the source term average, as follows:

S ± i = 1 ∆x x i+ 1 2 x i-1 2 S(p n w (x; i), x) dx + O(∆x δ ). (4.7)
In practice, this approximation is nothing but a quadrature formula of order δ, see for instance [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]. We then adopt the following well-balanced modification of this approximation:

S± i = 1 2 1 -θ n i-1 2 S n i-1 2 + 1 -θ n i+ 1 2 S n i+ 1 2 + 1 2 θ n i-1 2 + θ n i+ 1 2 S ± i . (4.8)
At this level, it is worth noting that the δ-order numerical scheme designed here is obtained arguing a very easy modification (4.6) and (4.8) of any polynomial reconstruction (4.2) of degree d and any source term integration (4.7) of order δ. However, this minor correction ensures that the scheme is well-balanced, is of order δ in space, and preserves the set of admissible states as soon as the original high-order scheme does.

Theorem 4. Given a well-balanced first-order scheme, the scheme (4.4), with the reconstructed states given by (4.5) and the source term approximation (4.8), satisfies the following properties:

(i) if w n
i is an approximation of w(x i , t n ) up to O(∆x δ ) for all i, then the scheme (4.4) is of order δ = d + 1 in space;

(ii) if (w n i ) i∈Z defines a steady state according to (2.3), then the scheme (4.4) is well-balanced, i.e. w n+1 i = w n i for all i ∈ Z;

(iii) the scheme (4.4) is robust, i.e. if the original reconstruction (4.2) of degree d preserves Ω, then Ω remains invariant by the well-balanced improvement of the reconstruction (4.6).

We stress that, according to Theorem 3, properties (i) and (ii) cannot be satisfied for perturbed steady states unless the initial condition is computed using (4.1) with, respectively, θ i = 1 for all i and θ i = 0 for all i.

Proof of Theorem 4. We establish properties (i), (ii) and (iii) in order.

(i) We first establish the order of accuracy, as defined in Theorem 1. To address such an issue, we consider w(x, t) a smooth unsteady solution of (1.1) so that the relation (3.9) holds. Next, with (w n i ) i∈Z given by (3.7) and (w n+1 i ) i∈Z given by the high-order scheme (4.4), the relation (3.8) holds for

F i+ 1 2 = f ± i+ 1 2 - 1 ∆t ∆t 0 f (w(x i+ 1 2 , t n + t)) dt, S i = S± i - 1 ∆t∆x ∆t 0 x i+ 1 2 x i-1 2 S(w(x, t n + t), x) dx dt,
where S± i is defined by (4.8) and f ± does not vanish, a Taylor expansion yields

θ n i+ 1 2 = 1 + O(∆x k ), with k ≥ δ = d + 1.
As a consequence, in the current unsteady context, by definition of the polynomial reconstruction according to (4.3), we obtain

w ± i+ 1 2 = w(x i+ 1 2 , t n ) + O(∆x δ ). Next, from (2.
2), we know that the numerical flux function is Lipschitz-continuous and consistent. Therefore, f ∆ (w -

i+ 1 2 , w + i+ 1 2 ) = f (w(x i+ 1 2 , t n )) + O(∆x δ
), and we get F i+ 1 2 = O(∆t) + O(∆x δ ). Concerning the order of accuracy of the source term, since (4.8) reduces to S± i = S ± i + O(∆x δ ), arguing (4.7) yields

S± i = 1 ∆x x i+ 1 2 x i-1 2 S(w(x, t n ), x) dx + O(∆x δ ). (4.9)
Plugging (4.9) into the definition of S i , we get

S i = O(∆t) + O(∆x δ ).
The establishment of the order of accuracy is thus completed.

(ii) For the proof of the well-balancedness property, let us consider (w n i ) i∈Z to define a steady state according to (2.3). By definition of the correction, given by (3.2) -(3.3), we easily obtain θ n i+ 1 2 = 0 for all i in Z so that the reconstructed states now read

w - i+ 1 2 = w n i and w + i+ 1 2 = w n i+1 .
Similarly, regarding the source term reconstruction given by (4.8), we now have S±

i = 1 2 (S n i-1 2 + S n i+ 1 2
). As a consequence, the high-order scheme (4.4) coincides with the first-order well-balanced scheme (2.1), and the preservation of the steady states immediately follows.

(iii) To conclude the proof, we now establish that the improvement (4.6) preserves the convex set Ω as long as the original polynomial reconstruction (4.2) preserves Ω. Indeed, we have

(p n w ) ± i = 1 -θ n i± 1 2 w n i + θ n i± 1 2 p n w x ± ∆x 2 ; i . Since p n w (x; i) ∈ Ω for all x ∈ (x i-1 2 , x i+ 1 2 ), w n i ∈ Ω and θ n i ∈ [0, 1], we immediately get (p n w ) ± i ∈ Ω.
The proof is thus completed.

Numerical experiments

To assess the performance of the scheme developed above, we now perform several numerical experiments. First, we describe in Section 5.1 the setup used to assess both the order of accuracy and the well-balancedness property of the schemes under consideration. Then, we apply the high-order well-balanced strategy to several systems, namely the shallow water equations with topography (1.4) in Section 5.2, the shallow water equations with friction (1.9) in Section 5.3, and the Euler equations with gravity (1.11) in Section 5.4. We also include, in Section 5.2.4, a discussion around the issues encountered with initialization.

Setup

To justify the relevance of the procedure outlined in Sections 3 and 4, we wish to compare the results of a given first-order well-balanced scheme to the ones produced by its second-order and high-order extensions, with and without the well-balancedness correction. For the sake of simplicity, we introduce the following notations:

• the P d scheme is the scheme of order d + 1 without the well-balancedness correction,

• the P WB d scheme is the scheme of order d + 1 with the well-balancedness correction.

Note that the P 0 and P WB 0 schemes are identical. Furthermore, note that forcing θ = 1 on the whole space-time domain in the P WB d scheme is enough to yield the P d scheme. In this paper, we consider high-order schemes up to a third-order accurate P 2 scheme. This is enough to justify both the high-order accuracy and the steady state preservation.

To use the P d scheme, we need to define three elements: the polynomial reconstruction from (4.2), the approximation of the source term average from (4.7), and the time integration. These elements are summarized in Table 5 means that the scheme is never well-balanced, while taking a very small C n i+ 1 2 means that the scheme is never high-order accurate. Investigating several expressions of C n i+ 1 2 led us to choose the numerical time derivative of the solution, since it is large when the solution is unsteady (and we want a large θ n i+ 1 2 ), and small when the solution is steady (and we want a small θ n i+ 1 2 ). As a consequence, it is the most generic choice we found that provided good numerical approximations for each system under consideration.

To that end, we define, for n ≥ 1,

C n i+ 1 2 = C θ 1 2 w n i+1 -w n-1 i+1 ∆t + w n i -w n-1 i ∆t ,
with C θ a constant parameter, which can be interpreted as a normalization of the time derivative. The choice of C θ depends on the numerical experiment under consideration (unless otherwise mentioned, we take C θ = 1). We also take C 0

i+ 1 2 = 1. Note that, equipped with this definition of C n i+ 1 2 , the expression (3.2) of θ n i+ 1 2 reads: θ n i+ 1 2 = ε n i+ 1 2 (C n i+ 1 2 ) k ε n i+ 1 2 (C n i+ 1 2 ) k + ∆x k . Therefore, we get θ n i+ 1 2 = 0 if ε n i+ 1 2 = 0 or if C n i+ 1 2
= 0. This is justified in each case, as follows.

• If ε n i+ 1 2 = 0, then a steady solution of the equations has been reached, and the first-order well-balanced scheme should be used to ensure the preservation of this solution. Taking θ n i+ 1 2 = 0 enables this behavior.

• If C n i+ 1 2 = 0, then a local steady solution of the P d scheme has been reached. Regardless of whether ε n i+ 1 2 = 0, we should get θ n i+ 1 2 = 0 in this case, since a steady solution for the P d scheme will not, in general, be a steady solution for the equations. Setting θ n i+ 1 2 = 0 for such cases perturbs the steady numerical solution and allows it to converge towards the true steady solution.

In practice, to avoid round-off errors, we take θ n i+ 1 2 = 0 as soon as the expression above yields a value smaller than 10 -12 .

This new expression of θ n i+ 1 2 gives a partial answer to the additional question of the interpretation of intermediate values of

θ n i+ 1 2 . Take θ n i+ 1 2 = ν, with 0 < ν < 1 some constant independent from ∆x. In this case, we get ε n i+ 1 2 (C n i+ 1 2 ) k = O(∆x k ). Since both quantities ε n i+ 1 2 and C n i+ 1 2
measure an error to a steady state, then at least one of these quantities is small, and the numerical solution is ∆x k -close to a steady state. Thus, heuristically, the first-order scheme is at least as accurate as the high-order scheme, since it is exact on steady states and the numerical solution is ∆x k -close to a steady state.

We shall perform three experiments for each system, in order to validate both the high-order accuracy and the well-balancedness property. These experiments are detailed below; system-specific parameters (such as the final physical time, for instance) will be given in the relevant sections.

In each experiment, the space domain is (0, 1) and the simulation is run until some final time t end . Each experiment relies on the following compactly supported C ∞ bump function:

ω(x) =                exp      1 - 1 1 -4 x - 1 2 2      if x - 1 2 < 1 4 , 0 otherwise.
In addition, unless otherwise mentioned, all errors computed in the remainder of the text are L2 errors between the approximate solution and the exact or reference solution.

The first experiment we perform yields a measure of the order of accuracy of the schemes, and it is designed to show that the well-balancedness correction does not reduce the accuracy for unsteady solutions. To correctly measure the order of accuracy, no slope limitation is added to the P d and P WB d schemes. Since we do not necessarily know an exact solution of the system under consideration, we compute a reference solution using a very fine grid made of 20 × 2 12 cells. Then, after computing the approximate solution on a coarser dyadic grid made of N = 20 × 2 k cells, 0 ≤ k 12, the fine solution is projected onto the coarser grid to measure the error between the reference solution and its approximation. To ensure that no shock waves form, the initial condition is smooth; its expression is given for each system. The initial condition is then evolved until the final time t end = 5 • 10 -3 . Periodic boundary conditions are prescribed for this experiment.

The second and third experiments assess the well-balancedness property. To that end, we study the dissipation of a perturbation applied to an initially steady solution. Here, we add a slope limitation to the P d and P WB d schemes, namely the MC limiter from [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] for d = 1 and the limiter from [START_REF] Schmidtmann | Relations Between WENO3 and Third-Order Limiting in Finite Volume Methods[END_REF] for d = 2. The initial condition is a steady solution w, computed by solving the nonlinear system (1.3). This steady solution is then perturbed using the bump function ω. Namely, each variable in w is multiplied by (1 + αω(x)). We distinguish two cases: a small and a large perturbation. For the small perturbation, we take α = β∆x 3 , with β a constant given for each system. For the large perturbation, we take α = 0.25. The steady solution is imposed on the boundaries, and the experiment is run until the numerical solution becomes steady; this final time is given for each system. We take 50 discretization cells for this experiment.

Application: the shallow water equations with topography

The first application concerns the shallow water equations with topography (1.4). The first-order wellbalanced scheme comes from [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF][START_REF] Berthon | A simple fully well-balanced and entropy preserving scheme for the shallow-water equations[END_REF]. It contains a parameter C, set here to +∞, or rather to the upper bound of the double-precision floating-point numbers in practice. In addition, the topography function is set to z(x) = ω(x) and we take g = 9.81.

Order of accuracy assessment

For this experiment, the initial condition is given by h 0 (x) = 2z(x) + cos 2 (2πx) and q 0 (x) = sin(2πx).

In Figure 5.1, we display the reference solution and the approximations given by the P 0 , P 2 and P WB 2 schemes with 40 discretization cells. We observe that the third-order schemes are very close to the reference solution, even with such a few cells.

This observation is confirmed in Table 5.2 and Figure 5.2, where we report the errors on h and q, as well as the orders of accuracy. As expected, the well-balancedness procedure does not alter the order of accuracy of the scheme, since the solution produced by the Table 5.4: Shallow water equations with topography: errors between the steady solution and the approximate solutions at time t end , for the steady state with a large perturbation, using 50 cells.

Well-balancedness property -small perturbation

The initial condition is computed the same way as in the previous test case, but using a small perturbation of amplitude β∆x 3 , with β = 10 -2 , instead of a large perturbation. We take C θ = 3.5 • 10 -2 , and the final physical time is t end = 5 • 10 -3 .

In Figure 5.5, we display the errors between the numerical solution and the underlying steady solution at the final time t end . The P 1 and P 2 schemes produce errors that destroy both the underlying steady state and the small perturbation. Indeed, these errors are actually larger than the perturbation, and the P 1 and P 2 schemes behave as though there were no perturbation. Conversely, the P 0 , P WB 1 and P WB 2 schemes yield comparable results, and the small perturbation is clearly visible against the underlying, exactly preserved steady state. Here, it was expected for the results of these three schemes to be close, since the convex combination procedure correctly switches from the locally less accurate high-order scheme to the locally more accurate well-balanced scheme. These observations are confirmed in Table 5.5, where the errors between the numerical solutions and the underlying steady states are reported. As expected, the well-balanced schemes produce an error of the order of magnitude of the perturbation amplitude (around 10 -7 ), while the errors made by using the non-well-balanced schemes are several orders of magnitude larger. We display the initial condition, as well as the approximations produced by the P 0 , P 2 and P WB 2 schemes at time t = 2.5 • 10 -2 , in Figure 5.9. We observe that the solutions of the P 2 and P WB 2 are indistinguishable, and that they are less diffusive than the solution given by the P 0 scheme. These observations are confirmed by the errors reported in Table 5 Table 5.9: Shallow water equations with Manning friction: errors between the underlying steady solution and the approximate solutions at time t end , for the steady state with a small perturbation, using 50 cells.

Application: the Euler equations with gravity

For the last application, we turn to another system to highlight the genericity of our method. We choose the Euler equations with gravity (1.11). The first-order well-balanced scheme is based on the strategy from [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography or Manning friction[END_REF][START_REF] Ghitti | A fully well-balanced scheme for the 1D blood flow equations with friction source term[END_REF]. It is designed to exactly preserve and capture the family of steady states given by (1.12).

We consider an ideal gas pressure law, where the pressure p is given by p

= (γ -1)(E - 1 2 q 2 ρ ).
Furthermore, we take the parameters γ = 1.4 and κ = 1. The gravity potential is given by φ(x) = ω(x).

Order of accuracy assessment

For this experiment, we take the following initial condition: ρ 0 (x) = 2 + cos 2 (2πx), q 0 (x) = sin(2πx), E 0 (x) = 5 + cos 2 (2πx).

In Figure 5.12, we display the reference solution and the approximations given by the P 0 , P 2 and P WB We also report the errors and the orders of accuracy for ρ, q and E in Table 5.10 and Figure 5.13. The same conclusion as for the shallow water system with topography or Manning friction is drawn. In Figure 5.15 and Table 5.12, the errors on ρ, q and E at the final time t end are reported, as well as the CPU time. The same conclusions as in the shallow water case are reached. We even observe a decrease in CPU time with the P WB 2 scheme compared to the P 2 scheme. This is due to the following interesting optimization. When a steady state is reached at interface x i+ 1 2 , we get θ i+ 1 2 = 0 and the high-order correction need not be computed, which saves some CPU time. x qq ex 0 0.5 Table 5.12: Euler equations with gravity: errors between the steady solution and the approximate at time t end , for the perturbed steady state experiment with 50 cells.

Well-balancedness property -small perturbation

This last experiment is a small perturbation of a steady solution, with parameter β = 10 -2 . We take C θ = 7 • 10 -2 , and the final physical time is t end = 5 • 10 -3 .

The errors between the approximate solution and the underlying steady state are displayed in Figure 5.16. Once again, like for the previous systems, we note that the P 1 and P 2 schemes have produced errors several orders of magnitude larger than the small perturbation, and destroyed both the perturbation and the underlying steady state. The numerical values of the errors reported in Table 5.5 once again validate this observation. x qq ex 0 0.5 Table 5.13: Euler equations with gravity: errors between the underlying steady solution and the approximate solutions at time t end , for the steady state with a small perturbation, using 50 cells.

Conclusion and perspectives

Usual high-order well-balanced methods involve costly inversions of nonlinear equations -or even systems of equations -to ensure both high-order accuracy and the preservation of complex steady solutions. The generic linear technique proposed in this manuscript avoids this extra cost while remaining able to exactly preserve steady states and to ensure high-order accuracy away from steady states.

Although this method has been presented in one space dimension, it can be naturally extended to more space dimensions. Indeed, consider the usual D-dimensional polynomial reconstruction of some variable w, which would replace the 1D polynomial reconstruction (4.2):

p n w (x; i) = w n i + π w i (x -x i ),
where x ∈ R D and π w i is a polynomial function of degree d, which satisfies some conditions detailed for instance in [START_REF] Diot | Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF][START_REF] Diot | The multidimensional optimal order detection method in the threedimensional case: very high-order finite volume method for hyperbolic systems[END_REF][START_REF] Michel-Dansac | A two-dimensional high-order well-balanced scheme for the shallow water equations with topography and Manning friction[END_REF]. To compute the high-order D-dimensional numerical flux, the above reconstruction is evaluated at several Gauss points on the cell interfaces. Let ζ i,j be a Gauss point at the interface between cell i and cell j. Then, in the spirit of (4.6), the evaluation of the polynomial at ζ i,j is modified as follows: p n w (ζ i,j ; i) := w n i + θ i,j π w i (ζ i,jx i ), where θ i,j is nothing but the error to the steady state between cells i and j. This D-dimensional extension is the object of future work.

Finally, numerical experiments show that the proposed technique performs as expected, i.e. we recover both high-order accuracy and steady state preservation at a negligible cost, for several hyperbolic systems.

The main advantages of this procedure thus are its genericity and its simplicity. Indeed, it is very easy to implement and, since it is linear, it does not involve much extra computational cost, if any. However, the technique requires prior knowledge of a first-order well-balanced scheme. Moreover, there remains one free parameter in the expression (3.2) of the steady state detector, but the value proposed in Section 5.1 seems relevant for each case treated so far.

An important difference between our method and existing literature, such as [START_REF] Castro | Well-balanced high order extensions of Godunov's method for semilinear balance laws[END_REF][START_REF] Castro Díaz | High order exactly well-balanced numerical methods for shallow water systems[END_REF][START_REF] Gómez-Bueno | High-order well-balanced methods for systems of balance laws: a control-based approach[END_REF], also needs to be mentioned. In these papers, a local steady solution is obtained in each cell by solving an ODE similar to (1.2); then, the deviations with respect to this local steady solution are reconstructed; and the final reconstruction is given by adding the reconstructed deviations to the local steady solutions. This reconstruction strategy, although much more computationally expensive than the one introduced in the present paper, is able to preserve the distance to a perturbation of a global steady solution. Indeed, a reconstruction of a perturbation of amplitude ε will remain ε-close to the underlying steady solution, which is not the case for our method. In our case, the perturbation is discretized with the high-order scheme, so with an error in O(∆x δ ), which may be much larger than ε for coarse meshes. Of course, this discussion becomes moot for fine meshes, since both strategies are consistent with a high order of accuracy.

( 3 . 9 )

 39 With (w n i ) i∈Z given by (3.7) and (w n+1 i ) i∈Z given by (3.5) -(3.6), a straightforward computation gives the expected relation (3.8), where we have set

i+ 1 2 2 ,

 22 and S±i are respectively given by (3.6) and (3.4). We first treat the approximation of the flux function. By definition of θ n i+1 given by (3.2), as long as ε n i+ 1 2

i+ 1 2 by ( 3 . 6 ) 2 , as long as ε n i+ 1 2

 23622 , with the high-order reconstructed states w ± establish the order of accuracy associated with the numerical flux function. First, because of the definition (3.2) of θ n i+1

2 Figure 5 . 2 :

 252 Figure 5.2: Shallow water equations with topography: error lines for the dyadic experiment. Left panel: error on h; right panel: error on q.

P0 scheme P1 scheme P WB 1 scheme

 1 scheme P2 scheme P WB 2 scheme error on h 4.28 • 10 -8 1.03 • 10 -3 4.27 • 10 -8 7.70 • 10 -4 4.54 • 10 -8 error on q 4.66 • 10 -8 7.94 • 10 -3 4.48 • 10 -8 6.63 • 10 -3 9.46 • 10 -8

2 Figure 5 . 8 :

 258 Figure 5.8: Shallow water equations with Manning friction: error lines for the dyadic experiment. Left panel: error on h; right panel: error on q.

2 schemes with 40 2 Figure 5 . 12 :

 2402512 schemes with 40 discretization cells. As usual, the third-order schemes are very close to the reference solution.

  2 ), and we get F i+1 2 = O(∆t) + O(∆x 2 ). Next, we study the accuracy of the source term discretization. By definition of the source term reconstruction (3.4), we have S± i = S ± i + O(∆x 2 ), and arguing (2.6) immediately yields S i = O(∆x 2 ) + O(∆t).

	Arguing Theorem 1, the second-order space accuracy is thus established.
	(ii) Concerning the preservation of the steady states, as soon as (w n i ) i∈Z satisfy (2.3), we easily get ε n i+ 1 2

Table 5 .

 5 1: Polynomial reconstruction from (4.2), source term average from (4.7), and time integrator for the P 1 and P 2 schemes.

	.1.
	Moreover, recall that the P WB d From this definition, one notes that larger values of C n scheme is defined up to the choice of C n i+ 1 2 i+ 1 2 increase θ n i+ 1 2 , while smaller values of C n in the definition (3.2) of θ n i+ 1 2 i+ 1 2 decrease θ n i+ 1 2 . Although Theorem 4 holds for any expression of C n , this choice will impact numerical . i+ 1 2 experiments. Indeed, if C n i+ 1 2 is too large, then θ n will also be too large, and steady states that should be i+ 1 2 captured will not be captured. Similarly, if C n i+ 1 2 is too small, θ n 2 i+ 1 will also be too small, and steady states

  Figure 5.1: Shallow water equations with topography: comparison between the reference and approximate solutions for the dyadic experiment with 40 discretization cells, at time t end = 5 • 10 -3 . Left panel: water height h; right panel: discharge q.

	h							q
								1
	2.5
								0
	2							-1
	0	0.2	0.4	0.6	0.8	1	x	0	0.2	0.4	0.6	0.8	1	x
					reference	P 0	P 2	P WB 2
		N	error, P0		order, P0	error, P2	order, P2 error, P WB 2	order, P WB 2
		40	1.04 • 10 -2	-	1.12 • 10 -3	-	1.12 • 10 -3	-
		80	5.24 • 10 -3	0.99	3.25 • 10 -4	1.78	3.26 • 10 -4	1.78
		160	2.58 • 10 -3	1.02	4.06 • 10 -5	3.00	4.08 • 10 -5	3.00
		320	1.29 • 10 -3	1.00	2.74 • 10 -6	3.89	2.73 • 10 -6	3.90
		640	6.42 • 10 -4	1.00	1.76 • 10 -7	3.96	1.76 • 10 -7	3.96
		1280 3.21 • 10 -4	1.00	1.34 • 10 -8	3.72	1.36 • 10 -8	3.69
		2560 1.60 • 10 -4	1.00	1.50 • 10 -9	3.16	1.54 • 10 -9	3.14
							P WB

Table 5 .

 5 2: Shallow water equations with topography: errors and orders of accuracy for the dyadic experiment. For the sake of conciseness, only the errors on h are reported in this table; the reader is referred to Figure5.2 for a visualization of the errors on q.

  Figure 5.4: Shallow water equations with topography: errors between the steady solution and the approximate solutions at time t end , for the steady state with a large perturbation, using 50 cells. Left panel: error on h; right panel: error on q.

	h -h ex				q -q ex
	10 -4
					10 -6
	10 -8
					10 -11
				x	10 -16	0	0.2	0.4	0.6	0.8	1	x
		P 0	P 1	P WB 1	P 2	P WB 2
		P0 scheme	P1 scheme P WB 1	scheme P2 scheme P WB 2	scheme
	error on h	5.72 • 10 -16 1.85 • 10 -3	1.08 • 10 -15	1.30 • 10 -3	1.91 • 10 -15
	error on q	2.15 • 10 -15 6.07 • 10 -3	2.64 • 10 -15	5.42 • 10 -3	4.73 • 10 -15
	CPU time (s)	2.91	8.59		9.5	23.78		24

Table 5 .

 5 5: Shallow water equations with topography: errors between the underlying steady solution and the approximate solutions at time t end , for the steady state with a small perturbation, using 50 cells. Figure 5.7: Shallow water equations with Manning friction: comparison between the reference and approximate solutions for the dyadic experiment with 40 discretization cells, at time t end = 5 • 10 -3 . Left panel: water height h; right panel: discharge q.

	h							q
	3
								1
	2.5							0
								-1
	0	0.2	0.4	0.6	0.8	1	x	0	0.2	0.4	0.6	0.8	1	x
					reference	P 0	P 2	P WB 2
		N	error, P0		order, P0	error, P2	order, P2 error, P WB 2	order, P WB 2
		40	1.04 • 10 -2	-	2.82 • 10 -4	-	2.82 • 10 -4	-
		80	5.23 • 10 -3	1.00	3.64 • 10 -5	2.95	3.64 • 10 -5	2.95
		160	2.57 • 10 -3	1.03	4.62 • 10 -6	2.98	4.62 • 10 -6	2.98
		320	1.28 • 10 -3	1.00	5.80 • 10 -7	2.99	5.80 • 10 -7	2.99
		640	6.40 • 10 -4	1.00	7.28 • 10 -8	3.00	7.28 • 10 -8	3.00
		1280 3.20 • 10 -4	1.00	9.11 • 10 -9	3.00	9.11 • 10 -9	3.00
		2560 1.60 • 10 -4	1.00	1.14 • 10 -9	3.00	1.14 • 10 -9	3.00

Table 5 .

 5 6: Shallow water equations with Manning friction: errors and orders of accuracy for the dyadic experiment. For the sake of conciseness, only the errors on h are reported in this table; the reader is referred to Figure5.8 for a visualization of the errors on q.

  Figure 5.9: Shallow water equations with Manning friction: comparison between the initial condition and the approximate solutions at time t = 2.5 • 10 -2 , for the perturbed steady state experiment with 50 cells. Left panel: water height h; right panel: discharge q.Figure 5.10: Shallow water equations with Manning friction: errors between the steady solution and the approximate solutions at time t end , for the perturbed steady state experiment with 50 cells. Left panel: error on h; right panel: error on q.Figure 5.11: Shallow water equations with Manning friction: errors between the underlying steady solution and the approximate solutions at time t end , for the steady state with a small perturbation, using 50 cells. Left panel: error on h; right panel: error on q.

	h -h ex			.7.					q -q ex
	10 -8	h									q 10 -6
	1.3										1.4 10 -9
	10 -11										1.2 10 -12
	1.2
	1.1 10 -14	0 0	0.2 0.2	0.4 0.4	0.6 0.6	P 0	0.8 0.8	P 1	1 1	x x P WB 1 10 -15 0 0.8 1	0	P 2	0.2 0.2	0.4 0.4 P WB 2	0.6 0.6	0.8 0.8	1	1 x	x
							initial	P 0		P 2		P WB 2
	h -h ex									q -q ex
	10 -7										10 -7
	10 -10										10 -10
	10 -13										10 -13
	10 -16	0	0.2	0.4	0.6		0.8		1	x	10 -16	0		0.2	0.4	0.6	0.8	1	x
						P 0		P 1		P WB 1			P 2		P WB 2
					P0 scheme	P1 scheme P WB 1	scheme P2 scheme P WB 2	scheme
			error on h 4.28 • 10 -11 1.54 • 10 -7		4.27 • 10 -11	3.09 • 10 -8	4.28 • 10 -11
			error on q 4.25 • 10 -11 3.60 • 10 -7		4.20 • 10 -11	7.70 • 10 -8	5.63 • 10 -11

Table 5 .

 5 10: Euler equations with gravity: errors and orders of accuracy for the dyadic experiment. For the sake of conciseness, only the errors on ρ are reported in this table; the reader is referred to Figure5.8 for a visualization of the errors on q and E. schemes at time t = 1 • 10 -1 , in Figure5.14. The solutions of the P 2 and P WB 2 schemes are once again indistinguishable and less diffusive than the one given by the P 0 scheme. As before, the errors reported in Table5.11 confirm this observation. Figure 5.14: Euler equations with gravity: comparison between the initial condition and the approximate solutions at time t = 1 • 10 -1 , for the perturbed steady state experiment with 50 cells. Left panel: density ρ; middle panel: momentum q; right panel: energy E.2.02 • 10 -2 1.15 • 10 -2 1.17 • 10 -2 6.65 • 10 -3 6.82 • 10 -3 error on q 3.34 • 10 -2 1.76 • 10 -2 1.60 • 10 -2 1.58 • 10 -2 1.56 • 10 -2 error on E 1.31 • 10 -1 7.06 • 10 -2 6.43 • 10 -2 5.13 • 10 -2 5.18 • 10 -2

	N	error, P0	order, P0	error, P2	order, P2 error, P WB 2	order, P WB 2
	80	4.18 • 10 -3	-	2.61 • 10 -4	-	2.61 • 10 -4	-
	160	2.07 • 10 -3	1.01	5.59 • 10 -5	2.22	5.59 • 10 -5	2.22
	320	1.03 • 10 -3	1.01	1.09 • 10 -5	2.36	1.09 • 10 -5	2.36
	640	5.14 • 10 -4	1.00	1.72 • 10 -6	2.66	1.72 • 10 -6	2.66
	1280 2.56 • 10 -4	1.00	2.31 • 10 -7	2.89	2.31 • 10 -7	2.89
	2560 1.28 • 10 -4	1.00	2.95 • 10 -8	2.97	2.95 • 10 -8	2.97
	5120 6.41 • 10 -5	1.00	3.70 • 10 -9	2.99	3.70 • 10 -9	2.99

Table 5 .

 5 11: Euler equations with gravity: L ∞ errors between the steady solution and the approximate solutions at time t = 1 • 10 -1 , for the perturbed steady state experiment with 50 cells.

  Figure 5.16: Euler equations with gravity: errors between the underlying steady solution and the approximate solutions at time t end , for the steady state with a small perturbation, using 50 cells. Left panel: error on ρ; middle panel: error on q; right error on E. .04 • 10 -8 2.38 • 10 -4 6.07 • 10 -8 1.00 • 10 -4 6.20 • 10 -8 error on q 4.90 • 10 -9 1.55 • 10 -3 4.69 • 10 -9 1.07 • 10 -3 2.79 • 10 -8 error on E 2.78 • 10 -7 1.83 • 10 -3 2.79 • 10 -7 7.13 • 10 -4 2.85 • 10 -7

					E -E ex
					10 -4
					10 -8
					10 -12
					10 -16	1	x
	P 0	P 1	P WB 1	P 2	P WB 2

P0 scheme P1 scheme P WB 1 scheme P2 scheme P WB 2 scheme error on ρ 6

scheme is almost the same as the one produced by the P 2 scheme in this unsteady context. We even observe a slight over-convergence, possibly explained by the use of the fourth-order accurate Simpson's method in the source term approximation.

Acknowledgments

C. Berthon acknowledges the support of ANR MUFFIN ANR-19-CE46-0004.

Well-balancedness property -large perturbation

The initial condition is a large perturbation of the steady solution implicitly given by (1.5), with

We take C θ = 9 • 10 -3 , and the final physical time is t end = 20.

In Figure 5.3, we display the initial condition, as well as the approximations given by the P 0 , P 2 and P WB 2 schemes at time t = 2.5 • 10 -2 . We observe that the solutions of the P 2 and P WB 2 are quite close, even in this case of a perturbed steady solution, and that they are less diffusive than the solution given by the P 0 scheme. To get a more precise error quantification, the errors between the approximate solutions and a reference solution computed with a fine mesh are presented in Table 5.3. We confirm the observations from Figure 5.3, since the well-balanced correction does not degrade the error produced by the corresponding high-order scheme. Note that small oscillations appear around the extrema with the third-order methods, both with and without the well-balanced correction. Such small oscillations are expected in this case, since the limiter from [START_REF] Schmidtmann | Relations Between WENO3 and Third-Order Limiting in Finite Volume Methods[END_REF] is less viscous than the MC limiter used in the second-order methods. Then, in Figure 5.4 and Table 5.4, we report the errors on h and q at the final time t end , as well as the CPU time taken by each method. We observe that the P 0 , P WB 1 and P WB 2 schemes have all converged towards the exact steady solution up to machine precision, while a non-zero error remains for the P 1 and P 2 schemes. These observations validate the proposed well-balancedness correction. We also note that, as expected, the well-balanced correction of the high-order methods is not computationally costly. Indeed, it corresponds to a 10% increase for the P WB 1 scheme and a 1% increase for the P WB 2 scheme. This can be contrasted to usual high-order well-balanced methods, where a costly nonlinear inversion is needed for the reconstruction. x qq ex P 0

.5: Shallow water equations with topography: errors between the underlying steady solution and the approximate solutions at time t end , for the steady state with a small perturbation, using 50 cells. Left panel: error on h; right panel: error on q.

Issues with initialization

For this last experiment on the shallow water system with topography (1.4), we wish to highlight the specific issue of initialization. discussed in Theorem 3. Indeed, let us underline that the scheme endowed with the correction is formally high-order accurate and exactly well-balanced. However, issues with the initialization procedure (4.1) mean that, on specific cases such as a perturbed steady state, the resulting discretization of the initial condition will either be exactly well-balanced but second-order accurate, or high-order accurate but non-well-balanced. We once again stress that the issues discussed in this paragraph stem from a failure of the initialization procedure, rather than a failure of the well-balanced correction itself.

For this experiment, we consider a large perturbation of a steady state. The initial condition is a perturbed steady solution, but this time the support of the perturbation is (3/8, 5/8) rather than (1/4, 3/4), which means that the steady solution is multiplied by (1 + αω(2x -1/2)). We compare the initialization given by (4.1) (called well-balanced initialization) to the high-order initialization, given by integrating the fine reference solution over each cell. To ensure that the solution remains smooth enough, we take t end = 10 -3 . In addition, we take C θ = 0.5.

The error lines are displayed in Figure 5.6. In the top left panel, the error on h clearly shows that the approximate solution obtained by the P 2 and P WB 2 schemes is capped to second order accuracy with the well-balanced initialization (4.1). The bottom left panel shows that using cell averages increases the order of accuracy to 3, as expected. Note that this behavior occurs whether or not the well-balanced correction is implemented: this confirms the fact that this is a failure of the initialization, rather than of the well-balanced correction. In the right panels, the same conclusions are drawn for q, although the effect is less visible. We also emphasize that, because the high-order scheme we consider in the experiments is third-order accurate, we only lose one order of accuracy by switching to the well-balanced initialization. For high-order schemes, more orders of accuracy would be lost on such test cases, and the discrepancy shown in Figure 5.6 would increase.

Application: the shallow water equations with Manning friction

The next application concerns the shallow water equations with Manning friction, governed by (1.9). The first-order well-balanced scheme comes from [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography or Manning friction[END_REF]. It also contains a parameter C, also set here to +∞. We take the friction exponent η = 7/3 according to Manning's model [START_REF] Manning | On the flow of water in open channels and pipes[END_REF][START_REF] Delestre | SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies[END_REF], and we set κ = 1 as well as g = 9.81.

N

L 2 error on q, HO initialization

.6: Shallow water equations with topography: error lines for the perturbed steady solution. Left panels: error on h; right panels: error on q. Top panels: well-balanced initialization with pointwise values; bottom panels: high-order accurate initialization with cell averages.

Order of accuracy assessment

The initial condition for this experiment is given by h 0 (x) = 2 + cos 2 (2πx) and q 0 (x) = sin(2πx).

In Figure 5.7, we display the reference solution and the approximations given by the P 0 , P 2 and P WB 2 schemes with 40 discretization cells. We once again observe that the third-order schemes are very close to the reference solution.

We report the errors on h and q, as well as the orders of accuracy, in Table 5.6 and Figure 5.8. The previous observation is once again confirmed since the P d and P WB d schemes produce almost exactly the same solution for this experiment.

Well-balancedness property -large perturbation

The initial condition is a perturbation of the steady solution implicitly given by (1.2), with

The final physical time is t end = 20.

P0 scheme P1 scheme P WB Table 5.7: Shallow water equations with Manning friction: L ∞ errors between the steady solution and the approximate solutions at time t = 2.5 • 10 -2 , for the perturbed steady state experiment with 50 cells.

In Figure 5.10 and Table 5.8, the errors on h and q at the final time t end are reported alongside the CPU time. Once again, the P 0 , P WB 1 and P WB 2 schemes have all converged towards the exact steady solution, and the P 1 and P 2 schemes produce a non-zero error. Like the previous test case, we note that the well-balanced correction is not costly.

Well-balancedness property -small perturbation

We consider a small perturbation of a steady solution, with β = 10 -5 ; the initial condition is computed in a similar fashion as in the large perturbation case. The final physical time is t end = 0.005, and we take C θ = 0.25.

The errors between the underlying steady solution and the numerical approximation are reported in Figure 5.11. As in the shallow water equations with topography, the large errors made by the P 1 and P 2 schemes destroy both the perturbation and the underlying steady state, as confirmed in Table 5.5.

Well-balancedness property -large perturbation

The initial condition is a perturbation of the steady solution implicitly given by (1.12), with

The final physical time is t end = 300. We display the initial condition, as well as the approximations produced by the P 0 , P 2 and P WB