Benjamin Legros
email: ben-jamin.legros@centraliens.net

Transient analysis of an affine Queue-Hawkes process

Keywords: Queueing, transient analysis, Queue-Hawkes, Laplace Transform

diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1 Introduction [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF] introduced the idea of self-excitement, a model in which the current intensity of events is determined by events in the past. In the so-called Hawkes processes, the rate of new event occurrences increases as each event occurs. The crowd-attraction behavior in Finance is well represented by Hawkes processes [START_REF] Errais | Affine point processes and portfolio credit risk[END_REF][START_REF] Fonseca | Hawkes process: Fast calibration, application to trade clustering, and diffusive limit[END_REF][START_REF] Haghighi | Modelling order arrivals at price limits using hawkes processes[END_REF]. Hawkes processes are also known to capture the spread of infectious diseases [START_REF] Kresin | Comparison of the Hawkes and SEIR models for the spread of Covid-19[END_REF], or the evolution of a population [START_REF] Brockwell | Birth, immigration and catastrophe processes[END_REF][START_REF] Kapodistria | Linear birth/immigration-death process with binomial catastrophes[END_REF][START_REF] Logachov | Large deviations in a population dynamics with catastrophes[END_REF]. In contact centers also, [START_REF] Altman | Do customer emotions affect worker speed? an empirical study of emotional load in online customer contact centers[END_REF] showed that the expressed sentiment of the customer influences agent response times and vice versa. The empirical analysis of [START_REF] Daw | The co-production of service: Modeling service times in contact centers using Hawkes processes[END_REF] further proved that this phenomenon of selfexcitement could be well modeled by Hawkes processes. In general, arrival processes with herding behavior in queueing systems can be modeled by Hawkes processes. For instance, [START_REF] Gao | Functional central limit theorems for stationary Hawkes processes and application to infinite-server queues[END_REF] analyzed a queue with Hawkes arrival process and infinite number of servers in the asymptotic regime where the baseline intensity is large.

The concept of Hawkes process has been extended to queueing processes via the definition of the Ephemerally Self-Exciting Process (ESEP) in [START_REF] Daw | The Queue-Hawkes process: Ephemeral self-excitement[END_REF]. In this paper, we analyze an ESEP with a linearly increasing arrival rate in the stationary and transient regimes. This model has been introduced in Section 3 of [START_REF] Daw | The Queue-Hawkes process: Ephemeral self-excitement[END_REF] and is called an affine Queue-Hawkes Process.

Specifically, we consider a single queue with a capacity of n customers and a pool of s homogeneous and independent servers, with 1 ≤ s ≤ n. Service times are exponentially distributed with rate 1. Customers arrive at the system according to a Poisson process with state-dependent arrival rate λ(x) = λ + γx, where

x is the number of customers present in the system, for 0 ≤ x < n. If a customer is not routed to service immediately upon arrival, then she/he waits in a queue for her/his turn to be served, with customers being served in order of arrival. One way to interpret this arrival process is to consider a primary flow of informed customers who are insensitive to the system state and a secondary flow of arriving customers generated by the presence of customers in the system. Each customer present in the system generates a new Poisson arrival process with rate γ. Note that our model can be adjusted to the evolution of a population with the parameters s and n being infinite, γ being a birth rate, 1 being a mortality rate, and λ being an immigration rate. For the transient analysis, we let N (t) be the number of customers in the system at time t ≥ 0 and set p x (t) = P (N (t) = x|N (0) = x 0 ), so that p x (t) is the transient probability to have x customers present in the system at time t given that the system is in state x 0 at t = 0, for 0 ≤ x 0 ≤ n. For the stationary regime, we denote by p x the probability to have x customers in the system. When n is finite, the number of states is finite. So, the system is stable. When n is infinite, the stability condition is s = ∞ and γ < 1 [START_REF] Daw | The Queue-Hawkes process: Ephemeral self-excitement[END_REF] or λ < s and γ = 0 as for an M/M/s queue.

Providing a transient analysis is useful for decision makers as the herding behavior seen in airports or in restaurants is often observed over short-time periods. Having a finite number of servers corresponds to situations where the system has a limited service capacity like in a restaurant or when a non-lethal disease can only be cured with hospital intervention. With an infinite number of servers, we mention [START_REF] Daw | Queues driven by Hawkes processes[END_REF] and [START_REF] Koops | Infinite-server queues with Hawkes input[END_REF] who determined the mean, and the moments of the Hawkes process-driven queue in both transient and steady state. Instead, in this paper, we express the Laplace Transform (LT) of the transient probabilities explicitly with finite or infinite number of servers and finite or infinite capacity. Having closed-form expressions is useful as it allows for a better understanding of the effect of the parameters and enables us to compute a specific probability without computing any other probabilities. Nevertheless, when the number of states is finite, other methods can be potentially more efficient to derive the transient probabilities. For instance, we mention the uniformization method introduced in [START_REF] Jensen | Markoff chains as an aid in the study of Markoff processes[END_REF] which was shown to be efficient in comparison with Runge-Kutta and Liou methods [START_REF] Grassmann | Transient solutions in Markovian queueing systems[END_REF] and in comparison with Padé's approximation [START_REF] Pulungan | Transient analysis of CTMCs: Uniformization or matrix exponential[END_REF].

For our queueing model, in states where some servers are idling, Laplace Transforms of transient probabilities are expressed using hypergeometric functions when γ ≥ 1. In the opposite case (as for the M/M/s queue), if the system starts when at least one server is idling, then the LT of transient probabilities are expressed using a contour integral. In states where all servers are busy, the LT of transient probabilities are expressed using the confluent hypergeometric functions of the first and second kind. Simplified expressions are provided when n or s tends to infinity. After applying a numerical inversion of the LT, we observe that the affine Queue-Hawkes process needs longer to reach stationarity compared to an M/M/s queue with the same flow of arrivals. This result is a consequence of over-dispersion of the arrival process. That is, the variance of the arrival process is larger than the mean, whereas the Poisson process has equal mean and variance. This phenomenon is also observed in other Hawkes processes [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF] and further justifies the need to obtain transient expressions.

In what follows, we review the literature related to the analysis of transient queues. The M/M/1 queue was the first transient queue analyzed [START_REF] Ledermann | Spectral theory for the differential equations of simple birth and death processes[END_REF][START_REF] Bailey | A continuous time treatment of a simple queue using generating functions[END_REF]. The transient queue-length distribution for this queue is expressed through modified Bessel functions of the first kind. Later, [START_REF] Abate | Transient behavior of the M/M/1 queue: Starting at the origin[END_REF] provided a transform factorization that helps developing approximations for the moments of the queue length in the M/M/1 queue. Several other approaches for the analysis of the M/M/1 queue have been developed. We refer to [START_REF] Van De Coevering | Computing transient performance measures for the M/M/1 queue[END_REF] for a review of the main findings for the computation of the performance measures of the M/M/1 queue. For the multi-server setting, [START_REF] Kelton | The transient behavior of the M/M/s queue, with implications for steady-state simulation[END_REF] investigated the transient behavior of the M/M/s queue and showed the implications of their results for simulation. As in this paper, they considered the issue of the speed of convergence to the stationary regime. Later, [START_REF] Parthasarathy | Transient solution to the many-server Poisson queue: a simple approach[END_REF] found a solution for the M/M/s queue from which the stationary regime could be derived.

Including abandonment or rejection renders the performance evaluation difficult. We mention [START_REF] Al-Seedy | Transient solution of the M/M/c queue with balking and reneging[END_REF] and [START_REF] Knessl | Transient analysis of the Erlang A model[END_REF] for the performance measures of the M/M/s+M queue, [START_REF] Ammar | The busy period of an M/M/1 queue with balking and reneging[END_REF] for the study of its busy period and [START_REF] Legros | Transient analysis of a Markovian queue with deterministic rejection[END_REF] when deterministic reneging is involved.

Structure of the paper. The rest of the paper is organized as follows. Section 2 presents the stationary behavior of the affine Queue-Hawkes process. Section 3 determines the Laplace Transform of the transient probabilities. Section 4 provides a numerical illustration for the computation of the transient probabilities.

The proofs are given in the appendix at the end of the paper.

Stationary analysis

We present the stationary behavior of our queueing model. We have (λ + γx)p x = min(x + 1, s)p x+1 , for 0 ≤ x < n. This leads to

p x = p 0 γ x Γ λ γ + x x!Γ λ γ , for 0 ≤ x ≤ s, and, p x = p 0 γ x Γ λ γ + x s!s x-s Γ λ γ , for s ≤ x ≤ n, with, (1)
p 0 =   s-1 x=0 γ x Γ λ γ + x x!Γ λ γ + n x=s γ x Γ λ γ + x s!s x-s Γ λ γ   -1 , (2)
where Γ(z) is the Gamma-function defined for all complex numbers except the non-positive integers as

Γ(z) = ∞ t=0 t z-1 e -t dt.
Recall that our queueing model has been studied in [START_REF] Daw | The Queue-Hawkes process: Ephemeral self-excitement[END_REF] with an infinite number of servers or with a finite number of servers and no queue (i.e., the blocking model). As such, the stationary probabilities in [START_REF] Daw | The Queue-Hawkes process: Ephemeral self-excitement[END_REF], given in Theorem 3.3 in the infinite number of servers case can be deduced from our results in Equations ( 1) and (2) by letting s tend to infinity under the stability condition γ < 1. Those given Table 1: Conditions for the monotonicity property of p x in x (s > 1 and n > s) Interval p x is increasing p x is decreasing p x has a minimum p x has a maximum

0 ≤ x ≤ s λ ≥ 1 and γ ≥ s-λ s-1 λ ≤ 1 and γ ≤ s-λ s-1 λ < 1 and γ > s-λ s-1 λ > 1 and γ < s-λ s-1 s ≤ x ≤ n γ ≥ s-λ s γ ≤ s-λ n-1 s-λ n-1 < γ < s-λ s -
in Section C for the blocking model can be obtained from Equations ( 1) and (2) by setting s = n. In Table 1, we specify the conditions on the system parameters to have p x increasing, decreasing, with one minimum, or with one maximum for the intervals 0 ≤ x ≤ s and s ≤ x ≤ n. This table is obtained by analyzing the sign of p x+1 px -1, for s > 1 and n > s. In the remaining case s = 1, we have p 1 ≥ p 0 if and only if λ ≥ 1. While [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF] did not consider a queueing model, they also show that the Hawkes process is over-dispersed.

This particular behavior of congestion-attraction queue, also viewed as an affine Queue-Hawkes process [START_REF] Daw | The Queue-Hawkes process: Ephemeral self-excitement[END_REF],

slows down the speed of convergence to the stationary regime (see Section 3).

Transient analysis

We dedicate this section to the transient analysis of this queue. We provide an analysis involving complex integral to represent the Laplace Transforms (LT) of the transient probabilities as in [START_REF] Knessl | Transient analysis of the Erlang A model[END_REF]. This method is an alternative to the usual method which consists of finding the LT of the transient probabilities from the computation of the probability-generating function (see for instance [START_REF] Gross | Fundamentals of Queueing Theory[END_REF] page 95) or the moment-generating function (see Section 3 of [START_REF] Daw | The Queue-Hawkes process: Ephemeral self-excitement[END_REF] for the affine Queue-Hawkes process similar to ours). This method enables us to directly express the solutions of the equations governing the evolution of the system state with known and tabulated confluent hypergeometric functions. The LT of the transient probabilities can be inverted either when the system starts with at least s customers or when the arrival rate generated by present customers (γ) is higher than or equal to the service rate.

The Markov chain associated with our queueing system can be split into two parts. In the first one, when the number of customers in the system is below s, the Markov chain behaves as in the evolution of a population with birth, death and immigration rates. The second part of the Markov chain, is the symmetric counterpart of the M/M/s+M queue. For this part of the Markov chain, the death rate is constant while the birth rate is linearly increasing. In the M/M/s+M queue, the birth rate is constant and the death rate is linearly increasing. Therefore, the method with contour integral adopted in [START_REF] Knessl | Transient analysis of the Erlang A model[END_REF] for the M/M/s+M queue provides us with interesting representations of the LT of the transient probabilities. As opposed to the M/M/s+M queue, the involved integrals are real and can be expressed with some known confluent hypergeometric functions.

Analysis of the Markov chain is provided in Section 3.1. Section 3.2 combines the results of the previous section to explicitly state the LT of the transient probabilities.

We start the analysis with the equations governing the evolution of the system state. The forward Kolmogorov equations are

p 0 (t) = p 1 (t) -λp 0 (t), (3)
p x (t) = (λ + (x -1)γ)p x-1 (t) -(λ + xγ + x)p x (t) + (x + 1)p x+1 (t), for 1 ≤ x ≤ s -1, p x (t) = (λ + (x -1)γ)p x-1 (t) -(λ + xγ + s)p x (t) + sp x+1 (t), for s ≤ x ≤ n -1, p n (t) = (λ + (n -1)γ)p n-1 (t) -sp n (t), with p x (0) = δ(x, x 0 ), where δ(x, x 0 ) = 1 if x = x 0 and δ(x, x 0 ) = 0, for x = x 0 .
We set the Laplace transform of the transient probabilities as

P * x (θ)=
∞ 0 e -θt p x (t) dt, the system of Equations ( 3) can be rewritten as

P * 1 (θ) -(θ + λ)P * 0 (θ) = -δ(0, x 0 ), (4)
(x + 1)P * x+1 (θ) + (λ + (x -1)γ)P * x-1 (θ) -(λ + x(γ + 1) + θ)P * x (θ) = -δ(x, x 0 ), for 1 ≤ x ≤ s -1, (5)
sP * x+1 (θ) + (λ + (x -1)γ)P * x-1 (θ) -(λ + xγ + s + θ)P * x (θ) = -δ(x, x 0 ), for s ≤ x ≤ n -1, (6)
(λ + (n -1)γ)P * n-1 (θ) -(s + θ)P * n (θ) = -δ(n, x 0 ). (7)
We next solve Equations ( 4)- [START_REF] Brockwell | Birth, immigration and catastrophe processes[END_REF], by considering the cases 0 < x 0 < s, s < x 0 < n, or x 0 = 0, s, n.

Remark: Equations ( 4)-( 7) can also be employed to obtain the LT of the first passage time from one state to another. From the LT of the first passage time, we may easily derive the moments of the passage time from successive derivatives.

3.1 Solution of Equations ( 5) and ( 6)

In Theorem 1, we give the solutions of Equations ( 5) and ( 6). These solutions are constructed from contour integral where a generating function is introduced. We determine the differential equation that this generating function should satisfy and find two independent solutions. We denote by F x (θ) and G x (θ) the two independent solutions of Equation ( 5) and by H x (θ) and I x (θ) those of Equation ( 6). The solutions of these equations are given as functions of the confluent hypergeometric function of the first and second kind, M (a, b, z) and U (a, b, z)

and the general hypergeometric function, Φ(a, b, c, z), for a, b, c, z ∈ C, Re(a) > 0, with |z| < 1 and c not equal to a negative integer, where

M (a, b, z) = ∞ k=0 Γ(a + k)Γ(b) Γ(a)Γ(b + k) z k k! , U (a, b, z) = 1 Γ(a) ∞ u=0 e -zu u a-1 (1 + u) b-a-1 du, and, Φ(a, b, c, z) = Γ(c) Γ(a)Γ(b) ∞ k=0 Γ(a + k)Γ(b + k) Γ(c + k)k! z k .
Theorem 1. We have

F x (θ) = x k=0 γ k k!(x -k)! Γ λ γ -θ 1-γ + k Γ λ γ -θ 1-γ Γ θ 1-γ + x -k Γ θ 1-γ
, for γ = 1,

F x (θ) = Γ(x + λ) Γ(x + 1)Γ(λ) M (-x, λ, -θ), for γ = 1, G x (θ) = γ -λ γ + θ γ-1 Γ x + λ γ Γ 1 + θ γ-1 Γ x + 1 + λ γ + θ γ-1 Φ λ γ + θ γ -1 , x + λ γ , x + 1 + λ γ + θ γ -1 , 1 γ , for γ > 1, G x (θ) = Γ(x + λ)U (x + λ, λ, θ), for γ = 1, G x (θ) = 1 2πi C z -(x+1) (z -1) -θ 1-γ (γz -1) -λ γ + θ 1-γ dz, for γ < 1, H x (θ) = Γ θ γ + 1 e -s γ U θ γ + 1, x + 1 + λ + θ γ , s γ , and
I x (θ) = M x + λ γ , x + 1 + λ + θ γ , - s γ Γ x + λ γ Γ 1 + θ γ Γ x + 1 + λ+θ γ ,
where C is a contour which starts at -∞ -i and ends at -∞ + i encircling z = 0, z = 1 γ , and z = 1 in the counterclockwise sense.

With this definition, we obtain F -1 (θ) = 0, F 0 (θ) = 1, and F 1 (θ) = θ + λ. Therefore, Equation ( 4) is also satisfied by F x (θ) but not by G x (θ). In the case of γ < 1, it is not possible to express G x (θ) explicitly or as a function of some known functions.

Model extension: λ(x) = λ + γ min(x, c), for 0 ≤ c ≤ n. The analysis with state-dependent arrival rate λ(x) = λ + γx considered in this paper can be extended to the case λ(x) = λ + γ min(x, c), for 0 ≤ c ≤ n to account for bounded arrival rates. The following intervals should be considered:

• For 0 ≤ x ≤ min(c, s), Equation ( 5) is unchanged. So, P x (θ) can be expressed as a linear combination of F x (θ) and G x (θ) as expressed in Theorem 1.

• If c ≤ s and c ≤ x ≤ s, then Equation ( 5) is modified by letting first γ tend to zero and next replace λ by λ + cγ. Observe that (1 -γt)

-λ γ + θ
1-γ tends to e λt as γ tends to zero. Therefore, from the analysis in the proof of Theorem 1, P x (θ) can be written as a linear combination of F x (θ) and G x (θ) defined as

F x (θ) = x k=0 (λ+cγ) k k!(x-k)! Γ(θ+x-k) Γ(θ) and G x (θ) = 1 2πi C z -(x+1) (z -1)
-θ e (λ+cγ)z dz.

• If c ≥ s and s ≤ x ≤ c, Equation ( 6) is unchanged. Thus, P x (θ) can be expressed as a linear combination of H x (θ) and I x (θ) as expressed in Theorem 1.

• For max(c, s) ≤ x ≤ n, the transitions are identical to those of an M/M/1 queue. Therefore, the results in [START_REF] Gross | Fundamentals of Queueing Theory[END_REF] page 99 apply and P x (θ) can be expressed as a linear combination of

Expressions of the transient probabilities

We are now in position to express the Laplace Transform of the transient probabilities. These expressions, given in Theorem 2 with λ(x) = λ + γx for 0 ≤ x ≤ n, involve the Wronskian of the sequences F x (θ), G x (θ), H x (θ)

and I x (θ). We define the Wronskian application W for two sequences U x (θ) and V x (θ), such that W

(U,V ) x (θ) = U x+1 (θ)V x (θ) -V x+1 (θ)U x (θ).
Theorem 2. We have with 0 ≤ x 0 < s,

P * x (θ) = W (G,H) s-1 (θ)-αn(θ)W (G,I) s-1 W (F,H) s-1 (θ)-αn(θ)W (F,I) s-1 (θ) F x 0 (θ) -G x 0 (θ) (x 0 + 1)W (G,F ) x 0 (θ) F x (θ), for, 0 ≤ x ≤ x 0 , P * x (θ) = F x 0 (θ) (x 0 + 1)W (G,F ) x 0 (θ) W (G,H) s-1 (θ) -α n (θ)W (G,I) s-1 (θ) W (F,H) s-1 (θ) -α n (θ)W (F,I) s-1 (θ) F x (θ) -G x (θ)
, for, x 0 ≤ x ≤ s, and,

P * x (θ) = F x 0 (θ) W (F,H) s-1 (θ) -α n (θ)W (F,I) s-1 (θ) γ s-x 0 -1 x 0 ! s! Γ λ γ + s -1 Γ λ γ + x 0 (H x (θ) -α n (θ)I x (θ)) , for, s ≤ x ≤ n.
For s ≤ x 0 ≤ n, we have

P * x (θ) = γ s s-x 0 -1 Γ λ γ + s -1 sΓ λ γ + x 0 α n (θ)I x 0 (θ) -H x 0 (θ) W (H,F ) s-1 (θ) -α n (θ)W (I,F ) s-1 (θ) F x (θ), for, 0 ≤ x ≤ s, P * x (θ) = α n (θ)I x 0 (θ) -H x 0 (θ) sW (I,H) x 0 (θ)(W (H,F ) s-1 (θ) -α n (θ)W (I,F ) s-1 (θ)) (W (I,F ) s-1 (θ)H x (θ) -W (H,F )
s-1 (θ)I x (θ)), for, s ≤ x ≤ x 0 , and,

P * x (θ) = I x 0 (θ)W (H,F ) s-1 (θ) -H x 0 (θ)W (I,F ) s-1 (θ) sW (I,H) x 0 (θ)(W (H,F ) s-1 (θ) -α n (θ)W (I,F ) s-1 (θ)) (H x (θ) -α n (θ)I x (θ)) , for, x 0 ≤ x ≤ n, where α n (θ) = sH n+1 (θ)-(λ+nγ)Hn(θ) sI n+1 (θ)-(λ+nγ)In(θ) .
In asymptotic cases, the formulas of Theorem 2 can be simplified. In what follows, we present the cases of n = ∞ and s = n = ∞.

Infinite capacity queue. We first consider the case with an infinite capacity queue (i.e., n = ∞). The analysis of Theorem 2 can be made in a similar way using lim x→∞

I x = 0 and lim x→∞ H x = ∞. With 0 ≤ x 0 ≤ s, we have P * x (θ) = F x 0 (θ)W (G,I) s-1 (θ) -G x 0 (θ)W (F,I) s-1 (θ) (x 0 + 1)W (G,F ) x 0 (θ)W (F,I) s-1 (θ) F x (θ), for 0 ≤ x ≤ x 0 , P * x (θ) = F x 0 (θ) (x 0 + 1)W (G,F ) x 0 (θ) W (G,I) s-1 (θ) W (F,I) s-1 (θ) F x (θ) -G x (θ)
, for x 0 ≤ x ≤ s, and,

P * x (θ) = F x 0 (θ) W (F,I) s-1 (θ) γ s-x 0 -1 x 0 !Γ λ γ + s -1 s!Γ λ γ + x 0 I x (θ), for x ≥ s.
With x 0 ≥ s, we get

P * x (θ) = I x 0 (θ) sW (F,I) s-1 (θ) γ s s-1-x 0 Γ λ γ + s -1 Γ λ γ + x 0 F x (θ), for 0 ≤ x ≤ s, P * x (θ) = I x 0 (θ) sW (H,I) x 0 (θ) H x (θ) - W (F,H) s-1 (θ) W (F,I)
s-1 (θ) I x (θ) , for s ≤ x ≤ x 0 , and,

P * x (θ) = W (F,I) s-1 (θ)H x 0 (θ) -W (F,H) s-1 (θ)I x 0 (θ) sW (F,I) s-1 (θ)W (H,I) x 0 (θ) I x (θ), for x ≥ x 0 .
Infinite number of servers. In the case with an infinite number of servers and an infinite capacity, the formulas can be further simplified into

P * x (θ) = G x 0 (θ) (x 0 + 1)W (F,G) x 0 (θ) F x (θ), for 0 ≤ x ≤ x 0 , and, P * x (θ) = F x 0 (θ) (x 0 + 1)W (F,G) x 0 (θ) G x (θ), for x ≥ x 0 .
In the case x 0 = 0 and γ = 1, we get

P * x (θ) = Gx(θ) (1+θ)G 0 (θ)-G 1 (θ) , where G x (θ) = ∞ t=0 e -θt t x-1+λ (1 + t) x+1 dt, for
x ≥ 0. Using an integration by part, we obtain

G 1 (θ) = -θ ∞ t=0 e -θt t λ 1 + t dt + λG 0 (θ). So, (λ + θ)G 0 (θ) -G 1 (θ) = θ ∞ t=0 e -θt t λ-1 dt = θ 1-λ Γ(λ). Finally, we obtain P * x (θ) = θ λ-1 Γ(λ) ∞ t=0 e -θt t x-1+λ
(1 + t) x+1 dt. In the case λ = 1, we deduce that p x (t) = t x (1+t) x+1 . This means that p x (t) has a geometric distribution. Since lim t→∞ p x (t) = 0, the combination of parameters λ = γ = 1 and s = ∞ leads to an unstable system.

Numerical illustration

The transient performance measures can be computed using a Laplace transform inversion. We use the speed up version of the Gaver-Stehfest algorithm presented in [START_REF] Cohen | Numerical methods for Laplace transform inversion[END_REF], page 144, equation (7.7), where a given function f (t) is approximated by ln(2)

t N j=1 K j •f * j ln(2) t
, where N is even and

K j = (-1) j+ N 2 min(j,N/2) k=[ j+1 2 ] k N/2 (2k)! (N/2-k)!k!(k-1)!(j-k)!(2k-j)! .
In Figures 2 and3, we derive p 0 (t) and p n (t) as functions of the time elapsed since the origin, t, for different values of the initial state. Figure 2 presents a situation where the stationary probabilities are strictly increasing in the system state while Figure 3 illustrates a case where extreme states have the highest stationary probabilities. When comparing Figures 2(a) and 2(b) with Figures 3(a) and 3(b), we observe that the speed of convergence to the stationary distribution is slower in Figure 3. In Figure 2, the state x = n = 15 plays a role of attraction in the evolution of the system state. In Figure 3, both states x = 0 and x = n = 15 play this role. This means in practice that the system state remains either around x = 0 or around x = 15 for a long time, which slows down the achievement of a close to stationary regime. This means that a system with congestion-attraction may take a long time before reaching the stationary regime as compared to classical queueing models used to represent service systems. It also means that system managers should be careful while using stationary performance measures when making staffing or routing decisions over a short-term horizon.

For instance, if the system initially starts empty, stationary performance measures could lead to over-staffing decisions if the congestion-attraction phenomenon is not observed for a sufficiently long time. Proof. We first consider Equation [START_REF] Larry | Special functions of mathematics for engineers[END_REF]. We propose to build two independent solutions to this equation. We introduce the sequence f x which satisfies (x + 1)

f x+1 + (λ + (x -1)γ)f x-1 = (λ + x(γ + 1) + θ)f x .
The idea is to propose to write f x as an integral and to determine the differential equation that the generating function F(z) should satisfy in order to generate some solutions to Equation [START_REF] Larry | Special functions of mathematics for engineers[END_REF]. From Cauchy' integral theorem, the function f x can be written as f x = C z -(x+1) F(z) dz, for some contour C and some function F(z). Then, we have

C z -(x+1) λ + x(γ + 1) + θ -(λ + (x -1)γ)z - x + 1 z F(z) dz = 0. (8)
We assume that C is such that there are no boundary contributions arising from endpoints of C in order to apply the integration by parts and write

xf x = C xz -(x+1) F(z) dz = C z -x F (z) dz. Therefore, Equation (8) leads to C z -(x+1) -(1 -z)(1 -γz)F (z) + (λ(1 -z) + θ) F(z) dz = 0. (9)
One solution of this equation can be obtained if for all z we have

(1 -z)(1 -γz)F (z) = (λ(1 -z) + θ) F(z). (10)
One solution of Equation ( 10) is

F(z) = (1 -z) -θ 1-γ (1 -γz) -λ γ + θ 1-γ , if γ = 1, and, F(z) = (1 -z) -λ e θ 1-z , if γ = 1.
Note that F(x) has branch points at z = 1 if γ < 1 and at z = 1/γ ifλ γ + θ 1-γ < 0. Based on F(z), we are in position to build two independent solutions for Equation (5) denoted by F x (θ) and G x (θ).

The function F x (θ).

Case 1: γ = 1. We consider the integral

F x (θ)= 1 2πi C 1 z -(x+1) (1 -z) -θ 1-γ (1 -γz) -λ γ + θ 1-γ dz, where C 1 is
a small circle in the z-plane on which |z| < min 1, 1 γ . The integrand in the expression of F x (θ) is analytic inside the circle C 1 . We can obtain an explicit expression of F x . Since (1 -z)

-θ 1-γ = ∞ k=0 Γ θ 1-γ +k Γ θ 1-γ k! z k , and (1 -γz)
-λ γ + θ 1-γ = ∞ k=0 Γ λ γ -θ 1-γ +k Γ λ γ -θ 1-γ k! γ k z k , we have (1 -z) -θ 1-γ (1 -γz) -λ γ + θ 1-γ = ∞ m=0 z m   m k=0 γ k k!(m -k)! Γ λ γ -θ 1-γ + k Γ λ γ -θ 1-γ Γ θ 1-γ + m -k Γ θ 1-γ   .
This leads to the expression of F x (θ) as in Theorem 1.

Case 2: γ = 1. We define F x (θ) as

F x (θ)= e -θ 2πi C 1 z -(x+1) (1 -z) -λ e θ 1-z dz,
where C 1 is defined as above.

The definition of F x (θ) corresponds to the integral representation of the Laguerre polynomials. Using Equation (10.50) page 399 in [START_REF] Larry | Special functions of mathematics for engineers[END_REF], we may then write

F x (θ) = Γ(x+λ) Γ(x+1)Γ(λ) M (-x, λ, -θ)
, where M (a, b, z) is the confluent hypergeometric function of the first kind.

The function G x (θ). Due to the initial condition at t = 0 we need another solution to Equation (5) (independent from F x (θ)). We differentiate the cases γ > 1, γ = 1, and γ < 1.

Case 1: γ > 1. In this case, we observe that F(1) = 0, we can then define G x as an integral on the real

interval (1, ∞) by G x (θ)= ∞ z=1 z -(x+1) (z -1) θ γ-1 (γz -1)
λ γ + θ γ-1 dz. We change the variable z by t = 1/z.

This leads to G x (θ)= 1 t=0 γ -λ γ + θ γ-1 t x-1+ λ γ (1 -t) θ γ-1 1 - t γ -λ γ + θ γ-1
dt. Therefore, G x can be related to the general hypergeometric function. This function has an integral representation when Re(c) > Re(b) > 0;

Φ(a, b, c, z)= Γ(c) Γ(b)Γ(c -b) 1 t=0 t b-1 (1 -t) c-b-1 (1 -zt) -a dt, which proves the expression of G x (θ).
Case 2: γ = 1. In this case, we observe that F(1 + ) = 0, we can then define G x as an integral on the real interval

(1, ∞) by G x (θ)= ∞ z=1 z -(x+1) (z -1) -λ e θ 1-z dz. We change the variable z in u = -1 1-z . This leads to G x (θ)= ∞ u=0 u x-1+λ (1 + u) -(x+1)
e -θu du. Thus, we may write G x (θ) = Γ(x + λ)U (x + λ, λ, θ).

Case 3: γ < 1. This case is more difficult since we may not find z such that F(z) = 0. We instead consider the complex integral with contour C which starts at -∞-i and ends at -∞+i encircling z = 0, z = 1 γ , and z = 1 in the counterclockwise sense. We thus define

G x as G x (θ)= 1 2πi C 2 z -(x+1) (z -1) -θ 1-γ (γz -1) -λ γ + θ 1-γ dz.
G x (θ) satisfies the homogeneous form of Equation ( 5).

We now consider Equation ( 6) with the same approach. We shall construct two independent solutions to this difference equation. Let g x satisfy sg x+1 + (λ + (x -1)γ)g x-1 = (λ + xγ + s + θ)g x and represent g x as a contour integral, with g x = C z -(x+1) G(z) dz, for some function G(.) and contour C. Then, we have

C z -(x+1) λ + xγ + s + θ -(λ + (x -1)γ)z - s z G(z) dz = 0. (11)
Again, we assume that C is such that there are no boundary contributions arising from endpoints of C in order to apply the integration by parts and obtain

C z -(x+1) γz(1 -z)G (z) + λ(1 -z) + s + θ - s z G(z) dz = 0.
One solution of the above equation can be obtained if for all z we have

γz(1 -z)G (z) = -λ(1 -z) + s + θ - s z G(z). (12)
One solution of Equation ( 12) is

G(z) = e -s γz z -λ+θ γ (1 -z) θ γ . (13)
The function G(z) has a unique branch point at z = 0. To determine two independent solutions of Equation ( 6), we shall find two intervals leading to two independent solutions. These functions are denoted by H x (θ)

and I x (θ).

The function H x (θ). Since G(1) = G(0 + ) = 0, we first consider the interval (0, 1). Thus, we define H x (θ)

as H x (θ) = 1 t=0 e -s γt t -(x+1)-λ+θ γ (1 -t) θ γ dt. With the change of variable u = 1 t -1, we obtain H x (θ) = e -s γ ∞ u=0 e -s γ u (u + 1)
x-1+ λ γ u θ γ du, which leads to the expression of H x (θ).

B Proof of Theorem 2

Proof. To simplify the writing, we do not express the dependence of the considered functions in θ in the proof.

We first consider the case 0 < x 0 < s. We introduce the constant A, B, C, D, and E such that P * x = AF x , for, 0 ≤ x ≤ x 0 , P * x = BF x + CG x , for, x 0 ≤ x ≤ s, and P * x = DH x + EI x , for, s ≤ x ≤ n.

The boundary equations allow us to relate the different constants. With Equation ( 7), we have (λ + (n - This leads to B = D (s+1)G s+1 (Hs-αnIs)-sGs(H s+1 -αnI s+1 ) (s+1)(FsG s+1 -GsF s+1 ) = Dβ 1 s,n , and C = -D (s+1)F s+1 (Hs-αnIs)-sFs(H s+1 -αnI s+1 ) (s+1)(FsG s+1 -GsF s+1 )

= -Dβ 2 s,n . We now use the condition x = x 0 . This leads to AF x 0 = D β 1 s,n F x 0 -β 2 s,n G x 0 , and, A(x 0 + 1)F x 0 +1 = D(x 0 + 1) β 1 s,n F x 0 +1 -β 2 s,n G x 0 +1 + 1.

Therefore, we have A = which further simplifies the expressions. The formulas for 0 < x 0 < s can be extended to the case x 0 = 0. The case s < x 0 < n can be obtained with a similar approach. We introduce the constant Ã, B, C, D, and Ẽ such that P * x = ÃF x , for, 0 ≤ x ≤ s, P * x = BH x + CI x , for, s ≤ x ≤ x 0 , and P * x = DH x + ẼI x , for, x 0 ≤ x ≤ n. One difference here is that the function G x is not involved in the solution. After using the boundary equations, the Wronskian of I x and H x is involved. We find that sW ( γ s )

x Γ λ γ +x Γ λ γ . The formulas obtained for s < x 0 < n are valid for x 0 = s and for x 0 = n.

Figure 1 Figure 1 :

 11 Figure1presents the stationary probabilities in two systems with 10 servers and a capacity of 15 customers. In the first system, we set λ = 1 and γ = 0.903 to represent a situation with a significant congestion-attraction, and in the second one we adopt λ = 6.803 and γ = 0 to represent a queue without the phenomenon of congestion-attraction. The parameters are chosen such that the stationary expected number of customers in the system, E(N ), is equal to 7 in both cases. In the queue without congestion-attraction, the stationary

Figure 2 :

 2 Figure 2: Numerical results (λ = 1, γ = 1.2, s = 10, n = 15, N = 50)

Figure 3 :

 3 Figure 3: Numerical results (λ = 0.4, γ = 1.1, s = 10, n = 15, N = 50)

1 u=0e

 1 The function I x (θ). We now define I x (θ) as I x (θ) = With the change of variable z = 1 t , we obtain I x (θ) = Since M (a, b, z) has the following integral representation M (a, b, z) = Γ(b) Γ(a)Γ(b -a) zu u a-1 (1 -u) b-a-1 du, for Re(b) > Re(a) > 0, we obtain the expression of I x (θ).

  1)γ)(DHn-1 + EI n-1 ) -(s + θ)(DH n + EI n ) = 0. Moreover, we have s(DH n+1 + EI n+1 ) + (λ + (n -1)γ)(DH n-1 + EI n-1 ) -(λ + nγ + s + θ)(DH n + EI n ) = 0.By subtracting these two equations, we gets(DH n+1 + EI n+1 ) -(λ + nγ)(DH n + EI n ) = 0. Therefore, E = -D sH n+1 -(λ + nγ)H n sI n+1 -(λ + nγ)I n = -Dα n .For x = s, we may write BF s + CG s = D(H s -α n I s ), and, (s + 1)(BF s+1 + CG s+1 ) = sD(H s+1 -α n I s+1 ).

  Fx 0 -β 2 s,n Gx 0 β 2 s,n (x 0 +1)(Fx 0 G x 0 +1 -Gx 0 F x 0 +1) , and, D =

Fx 0 β 2 s

 2 ,n (x 0 +1)(Fx 0 G x 0 +1 -Gx 0 F x 0 +1 ) . Using Equation (5), we find that (x+1)W


Vassilios A Dougalis

Dimitrios E Mitsotakis

Jean-Claude Saut

Jean-Claude Saut

Boussinesq

J.-C Saut

BOUSSINESQ SYSTEMS OF BONA-SMITH TYPE ON PLANE DOMAINS: THEORY AND NUMERICAL ANALYSIS

whether they are published or not. The documents may come

Introduction

In this paper we will study the Boussinesq system

η t + ∇ • v + ∇ • ηv -b∆η t = 0, v t + ∇η + 1 2 ∇|v| 2 + c∇∆η -b∆v t = 0, (1.1)
where b > 0 and c < 0 are constants. This system is the two-dimensional version of a system in one space variable originally derived and analyzed by Bona and Smith, [BS]. It belongs to a family of Boussinesq systems, [BCSI], [BCL], that approximate the three-dimensional Euler equations for the irrotational free surface flow of an ideal fluid over a horizontal bottom. In (1.1) the independent variables x = (x, y)

and t represent spatial position and elapsed time, respectively. The dependent variables η = η(x, t) and v = v(x, t) = (u(x, t), v(x, t)), represent quantities proportional, respectively, to the deviation of the free surface from its level of rest, and to the horizontal velocity of the fluid particles at some height above the bottom. In (1.1) the variables are nondimensional and unscaled and the horizontal velocity v is evaluated at a height z = -1 + θ(1 + η(x, t)), for some θ ∈ [ 2/3, 1]. (In these variables the bottom lies at height z = -1.) In terms of the parameter θ, the constants in (1.1) are given by the formulas b = (3θ 2 -1)/6 and c = (2 -3θ 2 )/3. (The system originally analyzed by Bona and Smith in [BS] corresponds to θ 2 = 1.)

The value θ 2 = 2/3 (i.e. c = 0) yields the BBM-BBM system, [BC], [DMSII], [Ch].

The system (1.1) is derived from the Euler equations, [BCSI], [BCL], under a long wavelength, small amplitude assumption. Specifically, one assumes that ε := A/h 0 ≪ 1, λ/h 0 ≫ 1 with the Stokes number S := Aλ 2 /h 3 0 being of O(1). (Here A is the maximum amplitude of the surface waves measured over the level of rest, z = -h 0 is the (constant) depth of the bottom, and λ is a typical wavelength.) If one takes S = 1, then, in nondimensional, scaled variables, appropriate asymptotic expansions in the Euler equations yield equations of the form

η t + ∇ • v + ε(∇ • ηv -b∆η t ) = O(ε 2 ), v t + ∇η + ε 1 2 ∇|v| 2 + c∇∆η -b∆v t = O(ε 2 ),
(1.2) from which (1.1) follows by unscaling to remove ε, and replacing the right-hand side by zero.

The Cauchy problem for (1.1) in the case of one spatial variable has been proved to be globally well posed for 2/3 < θ 2 ≤ 1 in appropriate classical and Sobolev space pairs, [BS], [BCSI]. The analogous problem for θ 2 = 2/3 is locally well posed, [BCSI]. In [DMSI] we considered a more general class of systems in the two-dimensional case and proved that the corresponding Cauchy problem is locally well posed in appropriate pairs of Sobolev spaces. Initial-boundary-value problems (ibvp's) for (1.1) for 2/3 ≤ θ 2 ≤ 1 on a finite interval in one space variable were analyzed in [ADMI], and in [BC] for θ 2 = 2/3. It was proved in [ADMI] that the ibvp with Dirichlet boundary conditions, wherein η and u are given functions of t at the endpoints of the interval, is locally well posed. The corresponding ibvp with reflection boundary conditions (η x = u = 0 at both endpoints of the interval) was shown in [ADMI] to be globally well posed; so is also the periodic ivp. In [DMSII] we analyzed three ibvp's for the BBM-BBM system on a smooth plane domain Ω, corresponding to homogeneous Dirichlet boundary conditions for η and v on ∂Ω, to homogeneous Neumann boundary conditions for η and v on ∂Ω, and to the (normal) reflective boundary conditions ∂η ∂n = 0, and v = 0 on ∂Ω, where n is the normal direction to the boundary. We showed that these ibvp's are well posed locally in time in the appropriate sense.

In Section 2 of the paper at hand we consider the Bona-Smith system (1.1) and pose it as an ibvp on a plane domain Ω under a variety of homogeneous boundary conditions on ∂Ω, including e.g. homogeneous Dirichlet b.c.'s for η and v and reflective b.c.'s. We prove that the corresponding ibvp's are well posed, locally in time.

Turning now to the numerical solution of ibvp's for systems of the type (1.1) by Galerkin-finite element methods, we note first that it is quite straightforward to construct and analyze such schemes for the BBM-BBM system. For example, in [DMSI] we proved optimal-order L 2 -error estimates for the standared Galerkin semidiscretization of the BBM-BBM system with homogeneous boundary conditions on a smooth domain with a general triangulation. When c > 0, i.e. in the case of the proper Bona-Smith systems, the presence of the term ∇∆η complicates issues. In [DMSI] we analyzed the standard Galerkin semidiscretization with bicubic splines for this class of systems posed on rectangles with homogeneous Dirichlet boundary conditions, and proved optimal-order H 2 -error estimates for the approximation of η.

(Experimental evidence indicates that the L 2 -errors for the approximation of η, u and v with this scheme have suboptimal -O(h 3 ) -rate of convergence. In the one-dimensional case one may derive optimal-order estimates for the approximations of

η, u, v in W 1,∞ × L ∞ × L ∞ , cf. [ADMII].)
In Section 3 of the present paper we consider the Bona-Smith systems with c 0 posed on convex smooth planar domains with reflective boundary conditions. The systems are discretized on an arbitrary triangulation of the domain using a modified Galrkin method, wherein the Laplacian in the ∇∆η terms in (1.1) is discretized weakly by an appropriate discrete Laplacian operator. This enables us to prove optimal-order L 2 -and H 1 -error estimates on finite element subspaces of continuous, piecewise polynomial functions that include the case of piecewise linear functions utilized in most applications. The systems of ordinary differential equations representing the semidiscretizations of the Bona-Smith systems are shown to be non-stiff. One may thus use any explicit method for their temporal discretization.

We close the paper by showing the results of a series of numerical experiments of simulations of surface waves in complex domains, aimed at comparing the numerical solution of a Bona-Smith system with analogous results obtained by solving the BBM-BBM system.

2 Well-posedness of ibvp's for the Bona-Smith system Let Ω be a bounded plane domain with smooth boundary (or a convex polygon). We consider the Bona-Smith system

η t + ∇ • (v + ηv) -b∆η t = 0,
(2.1a) the norm and inner product, respectively, of L 2 = L 2 (Ω), by • s the norm of H s , and by

v t + ∇ η + 1 2 |v| 2 + c∆η -b∆v t = 0, (2.1b) for (x, t) ∈ Ω × R + , with b > 0, c < 0, with initial conditions η(x, 0) = η 0 (x), v(x, 0) = v 0 (x) x ∈ Ω. ( 2
• ∞ the norm of L ∞ = L ∞ (Ω).
The boundary conditions on η will be of the form

Bη = 0, x ∈ ∂Ω, t ∈ R + , (2.3a)
where the linear operator B and the domain Ω will be assumed to be such that the boundary-value problem

     -b∆w + w = f, in Ω, Bw = 0, in ∂Ω,
has for each f in L 2 a unique solution w ∈ H 2 for which Bw| ∂Ω is well defined. (We will also assume that an H 1 solution w of the problem is defined whenever f ∈ H -1 ). The boundary conditions on v will be of homogeneous Dirichlet type, i.e. = 0, η| Γ2 = 0.

v = 0, x ∈ ∂Ω, t ∈ R + . ( 2
In the sequel we let X := H 2 (Ω) ∩ {w : Bw = 0 on ∂Ω} and H 1 0 = (H 1 0 ) 2 etc. The main result of this section is Theorem 2.1 Given η 0 ∈ X, v 0 ∈ H 1 0 , there exists T > 0 and a unique solution (η, v) ∈ C([0, T ]; X) ∩ C([0, T ]; H 1 0 ) of the ibvp (2.1), (2.2), (2.3). Moreover, for each integer k 0

∂ k η ∂t k , ∂ k v ∂t k ∈ C([0, T ]; X) ∩ C([0, t]; H 1 0 ).
Proof: Write (2.1a) and (2.1b) as

η t + (I -b∆) -1 ∇ • (v + ηv) = 0, (2.4a) v t + (I -b∆) -1 ∇ η + 1 2 |v| 2 + c∆η = 0, (2.4b)
for x ∈ Ω, t > 0. In (2.4a) (I -b∆) -1 denotes the inverse of the operator I -b∆ with domain X, while in (2.4b) (I -b∆) -1 represents the inverse of I -b∆ with domain H 2 ∩ H 1 0 . Let F be the vector field on X × H 1 0 defined by

F (η, v) := (I -b∆) -1 ∇ • (v + ηv), (I -b∆) -1 ∇(η + 1 2 |v| 2 + c∆η) .
F is well defined on X × H 1 0 , since, by the Sobolev imbedding theorem, ηv ∈ H 1 and |v| 2 ∈ L 2 . Hence,

(I -b∆) -1 ∇ • (v + ηv) ∈ X and (I -b∆) -1 ∇(η + 1 2 |v| 2 + c∆η) ∈ H 1 0 . Moreover, F is C 1 on X × H 1 0 ,
with derivative F ′ (η * , v * ) given by

F ′ (η * , v * )(η, v) = (I -b∆) -1 ∇ • (v + ηv * + η * v), (I -b∆) -1 ∇(η + v * • v + c∆η) .
The continuity of F ′ follows from the Sobolev imbedding theorem and the regularity properties of the operators (I -b∆) -1 , considered as inverses of (I -b∆) on X for the first component, and of I -b∆ on

H 2 ∩ H 1 0 for the second.
By the standard theory of ordinary differential equations in Banach spaces, we conclude therefore that there exists a unique maximal solution (η, v)

∈ C 1 ([0, T ]; X) × C 1 ([0, T ]; H 1 0 ) of (2.4) for some T > 0,
with η| t=0 = η 0 and v| t=0 = v 0 . The first conclusion of the theorem follows. The assertion on

∂ k η ∂t k , ∂ k v ∂t k
follows by differentiating (2.4) with respect to t k -1 times.

Remark 2.1 It is not hard to see, using the energy method on the nondimensional but scaled system

(1.2) (with right-hand side replaced by zero), that at least in the cases of homogeneous Dirichlet boundary conditions for η and v and reflective boundary conditions ∂η ∂n = 0, v = 0 on ∂Ω, the maximum existence time T ε is independent of ε.

Remark 2.2 If the domain Ω is smooth, one gets smooth solutions from smooth data. Namely, let k ∈ N, k ≥ 3, and assume that

η 0 ∈ X ∩ H k , v 0 ∈ H 1 0 ∩ H k-1 , then (η, v) ∈ C([0, T ]; X ∩ H k ) ∩ C([0, T ]; H 1 0 ∩ H k-1
). This follows directly from the elliptic regularity estimates on I -b∆ (with the ad hoc boundary conditions). One is thus reduced to an ODE in

X ∩ H k × H 1 0 ∩ H k-1 . Remark 2.3 Consider the ibvp (2.1)-(2.
2) for the Bona-Smith system with (normal) reflective bound-

ary conditions ∂η ∂n = 0, v = 0, for x ∈ ∂Ω, t ∈ R + . (2.5)
It is known that in one dimension, cf. [ADMI], the Hamiltonian of the ibvp (2.1), (2.2), (2.5),

H = H(η, v) = 1 2 Ω η 2 + (1 + η)|v| 2 -c|∇η| 2 dx, (2.6)
is conserved. (Indeed this implies global existence-uniqueness of the solution of this ibvp in 1D provided that η 0 (x) + 1 > 0 and H(η 0 , u 0 ) is suitably restricted.) In the two-dimensional case H is not conserved in general. To see this, write the system (2.1) for (x, t) ∈ Ω × R + as

η t + ∇ • P = 0, (2.7a) v t + ∇Q + b∇ × ω t = 0, (2.7b)
where

P := v + ηv -b∇η t , Q := η + 1 2 |v| 2 + c∆η -b∇ • v t ,
and ω is the vorticity of the flow given by

ω = ∇ × v = (0, 0, ω), ω := v x -u y . From (2.7) we obtain Ω [η t Q + v t • P + ∇ • (QP) + bP • (∇ × ω t )] dx = 0 (2.8)
Using now the reflective boundary conditions in (2.8) and integrating by parts we see that in the maximal temporal interval of existence of a solution of the ibvp (2.1), (2.2), (2.5),

dH dt + b Ω P • (∇ × ω t )dx = 0. (2.9)
Hence, a simple sufficient condition for the conservation of the Hamiltonian is irrotationality of the flow.

In one space dimension, the flow is trivially irrotational. In 2D taking the curl of (2.1a) we see that

∂ t (ω -b∆ω) = 0, t > 0.
(2.10) Now, if the flow is, for example, irrotational at t = 0, the reflective boundary conditions do not allow us to conclude from (2.10) that ω = 0 for t > 0. However, in the case of the Cauchy problem or the ibvp with periodic boundary conditions on η and v on a rectangle (wherein (2.8) holds as well), ω(0) = 0 implies by (2.10) that ω = 0 for t > 0. (This was also noticed for the BBM-BBM system in [CI].) In these cases, it follows by (2.9) that the Hamiltonian is invariant. Note however that conservation of H does not imply global well-posedness in 2D.

A modified Galerkin method

We turn now to the numerical solution of the ibvp (2.1), (2.2), (2.5) by Galerkin-finite element methods.

The usual Galerkin method for Bona-Smith systems that was analyzed in [DMSI] requires finite element subspaces consisting of C 2 functions. Hence, it is not very useful in practice, as it cannot be applied to arbitrary plane domains and triangulations. In this section a modified Galerkin method for the numerical solution of the ibvp with (normal) reflection boundary conditions is analyzed. The method is more versatile, being applicable to general triangulations and domains (even with piecewise linear continuous functions) and is shown to have error estimates with optimal convergence rates in L 2 and H 1 .

For ease of reference we rewrite here the ibvp that we will approximate. We seek η and v, defined for (x, t) ∈ Ω × [0, T ] and satisfying for constants b > 0 and c ≤ 0

η t + ∇ • v + ∇ • ηv -b∆η t = 0, v t + ∇(η + 1 2 |v| 2 + c∆η) -b∆v t = 0, (x, t) ∈ Ω × [0, T ], η(x, 0) = η 0 (x), v(x, 0) = v 0 (x), x ∈ Ω, ∂η ∂n (x, t) = 0, v(x, t) = 0, (x, t) ∈ ∂Ω × [0, T ]. (R)
We assume that Ω is a convex domain with smooth enough boundary ∂Ω and that the ibvp (R) has a unique solution (η, v) which is smooth enough for the purposes of its numerical approximation. In the sequel, we put x = (x, y).

We suppose that T h is a regular triangulation of Ω with triangles τ of maximum sidelength h and let S h be a finite-dimensional subspace of C( Ω) ∩ H 1 , which, for small enough h and integer r ≥ 2, satisfies the approximation property inf

χ∈ e S h { w -χ + h w -χ 1 } ≤ Ch s w s , 1 ≤ s ≤ r, (3.1)
when w ∈ H s . (In this section C will denote generic constants independent of h.) On the same triangulation we denote by S h the subspace of S h consisting of the elements of S h that vanish on the boundary i.e. S h = S h ∩ H 1 0 . Hence, on S h we have inf

χ∈S h { w -χ + h w -χ 1 } ≤ Ch s w s , 1 ≤ s ≤ r, (3.2)
for w ∈ H s ∩ H 1 0 . We will assume that the elements of S h and S h are piecewise polynomial functions defined on T h , of degree at most r -1 on each τ ∈ T h .

We consider the symmetric bilinear form a D :

H 1 0 × H 1 0 → R defined by a D (u, v) := (u, v) + b(∇u, ∇v), u, v ∈ H 1 0 , (3.3) which is coercive on H 1 0 × H 1 0 .
In addition, we consider the symmetric bilinear form a N :

H 1 × H 1 → R that is defined by a N (u, v) := (u, v) + b(∇u, ∇v), u, v ∈ H 1 , (3.4)
and is coercive on

H 1 × H 1 .
With the aid of a D , a N we define the elliptic projection operators R h :

H 1 0 → S h , R h : H 1 → S h as follows: a D (R h w, χ) = a D (w, χ), ∀w ∈ H 1 0 , χ ∈ S h , (3.5a) a N ( R h w, χ) = a N (w, χ), ∀w ∈ H 1 , χ ∈ S h . (3.5b)
As a consequence of (3.1), (3.2) and elliptic regularity we have then

w -R h w k ≤ Ch s-k w s , w ∈ H s ∩ H 1 0 , 2 ≤ s ≤ r, k = 0, 1, (3.6a) and w -R h w k ≤ Ch s-k w s , w ∈ H s , 2 ≤ s ≤ r, k = 0, 1. (3.6b)
We assume that the triangulation T h is quasiuniform. Then, the following inverse assumptions hold on

S h , [Ci] χ 1 ≤ Ch -1 χ , ∀χ ∈ S h , (3.7) χ ∞ ≤ Ch -1 χ , ∀χ ∈ S h , (3.8)
The same inverse assumptions hold on S h as well.

We will also assume that for the elliptic projections we have the following approximation properties in the L ∞ norm:

w -R h w ∞ ≤ Cγ(h) w r,∞ , ∀w ∈ W r ∞ ∩ H 1 0 , (3.9a) w -R h w ∞ ≤ Cγ(h) w r,∞ , ∀w ∈ W r ∞ , (3.9b)
where γ(h) = h r | log h| r with r = 0 if r > 2 and r = 1 when r = 2, cf. [S]. Here, W r ∞ (with norm • r,∞ )

denotes the L ∞ -based Sobolev space on Ω of order r.

In the sequel we will use the discrete Laplacian operator ∆ h : H 1 → S h , defined for w ∈ H 1 by ( ∆ h w, χ) = -(∇w, ∇χ), ∀χ ∈ S h .

(3.10)

By the divergence theorem, it is easy to check that for w ∈ H 2 there holds

( ∆ h R h w, χ) = (∆w, χ) - ∂Ω ∂w ∂n χds - 1 b (w, χ) + 1 b ( R h w, χ), χ ∈ S h . (3.11) Hence, if w ∈ H 2 with ∂w ∂n | ∂Ω = 0, then ∆ h R h w = P h ∆w + 1 b ( R h -P h )w, (3.12)
where P h is the L 2 -projection onto S h . The identity (3.11) is proved as follows: For w ∈ H 2 , χ ∈ S h , using the definition of R h we have

( ∆ h R h w, χ) = -(∇ R h w, ∇χ) - 1 b ( R h w, χ) + 1 b ( R h w, χ) = -(∇w, ∇χ) - 1 b (w, χ) + 1 b ( R h w, χ), = - ∂Ω ∂w ∂n χds + (∆w, χ) - 1 b (w, χ) + 1 b ( R h w, χ).
Denoting the components of v as (u, v), we define the semidiscrete modified Galerkin method as follows.

We seek

η h : [0, T ] → S h , u h , v h : [0, T ] → S h , approximations to η, u, v, respectively, such that a N (η ht , φ) + (u hx , φ) + (v hy , φ) + ((η h u h ) x , φ) + ((η h v h ) y , φ) = 0, φ ∈ S h , a D (u ht , χ) + (η hx , χ) + (u h u hx , χ) + (v h v hx , χ) + c(( ∆ h η h ) x , χ) = 0, χ ∈ S h , a D (v ht , ψ) + (η h y , ψ) + (u h u hy , ψ) + (v h v hy , ψ) + c(( ∆ h η h ) y , ψ) = 0, ψ ∈ S h , η h (•, 0) = R h η 0 , u h (•, 0) = R h u 0 , v h (•, 0) = R h v 0 . 0 ≤ t ≤ T, (3.13)
These relations are discrete analogs to the corresponding variational forms of the first p.d.e. of (R) in

H 1 and of the second in (H 1 0 ) 2 . Note that since η h u h ∈ C( Ω), η h u h | τ ∈ C ∞ (τ ) for each τ ∈ T h , and
η h u h | ∂Ω = 0, it follows that η h u h ∈ H 1 0 . Similarly, η h v h ∈ H 1 0 , u 2 h ∈ H 1 0 , v 2 h ∈ H 1 0 .
Hence, all terms in the inner products of (3.13) are well defined. We now consider the mappings fx , fy : L 2 → S h defined for w ∈ L 2 by:

a D ( fx (w), φ) = (w, φ x ), φ ∈ S h , a D ( fy (w), φ) = (w, φ y ), φ ∈ S h ,
and the mappings and ĝx .ĝ y :

L 2 → S h defined for w ∈ L 2 by a N (ĝ x (w), χ) = (w, χ x ), χ ∈ S h , a N (ĝ y (w), χ) = (w, χ y ), χ ∈ S h .
Then (3.13) can be written as

η ht = F (u h , v h , η h ), u ht = G(u h , v h , η h ), 0 ≤ t ≤ T v ht = Z(u h , v h , η h ), η h (0) = R h η 0 , u h (0) = R h u 0 , v h (0) = R h v 0 , (3.14)
where

F (u h , v h , η h ) := ĝx (u h ) + ĝy (v h ) + ĝx (η h u h ) + ĝy (η h v h ), G(u h , v h , η h ) := fx (η h ) + 1 2 ( fx (u 2 h ) + fx (v 2 h )) + c fx ( ∆ h η h ), Z(u h , v h , η h ) := fy (η h ) + 1 2 ( fy (u 2 h ) + fy (v 2 h )) + c fy ( ∆ h η h ).
For the mappings fx , fy , ĝx and ĝy we have the following stability estimates:

Lemma 3.1 There exists a constant C such that

(i) fx (w) 1 ≤ C w , and fy (w) 1 ≤ C w , w ∈ L 2 .
(ii) ĝx (w) 1 ≤ C w , and ĝy (w

) 1 ≤ C w , w ∈ L 2 .
Proof: The proof follows immediately from the coercivity of a D on H 1 0 × H 1 0 and of a N on H 1 × H 1 and the definitions of fx , fy , ĝx , ĝy .

We define now the negative norms • -1 and • -2 for functions w ∈ L 2 as w -1 := sup

z∈H 1 z =0
(w, z) z 1 and w -2 := sup

z∈H 2 ∩H 1 0 z =0 (w, z) z 2 .
Lemma 3.2 There exists a constant C > 0 such that

fx (χ) ≤ C χ -1 , fy (χ) ≤ C χ -1 , ∀χ ∈ S h . (3.15)
Proof: Let χ ∈ S h . Consider the problem Lw = χ x with w = 0 on ∂Ω, where L := I -b∆ with domain

H 2 ∩ H 1 0 .
Then, by elliptic regularity and denoting by L -1 the inverse of L, we have

χ x -2 = sup 0 =z∈H 2 ∩H 1 0 (χ x , z) z 2 = sup 0 =z∈H 2 ∩H 1 0 (Lw, z) z 2 ≥ (Lw, L -1 w) L -1 w 2 = w 2 L -1 w 2 ≥ C w 2 w = C w . (3.16) Let now 0 = z ∈ H 2 ∩ H 1 0 . Then (χ x , z) z 2 = -(χ, z x ) z 2 ≤ χ -1 z x 1 z 2 ≤ χ -1 . Hence sup 0 =z∈H 2 ∩H 1 0 (χ x , z) z 2 ≤ χ -1 . Therefore χ x -2 ≤ χ -1 ,
and by (3.15)

w ≤ C χ -1 .
(3.17)

Consider R h w, the elliptic projection of w onto S h . Note that R h w = fx (χ). In addition, by (3.6a)

R h w -w ≤ Ch 2 w 2 ≤ Ch 2 χ x ≤ Ch 2 h -2 χ -1 = C χ -1 .
(In the last inequality we used the inverse inequality ϕ 1 ≤ Ch -2 ϕ -1 , ∀ϕ ∈ S h , which is valid since by (3.7):

ϕ -1 = sup 0 =z∈H 1 (ϕ, z) z 1 ≥ ϕ 2 ϕ 1 ≥ Ch 2 ϕ 2 1 ϕ 1 = Ch 2 ϕ 1 .)
We conclude that

R h w -w ≤ R h w -w ≤ C χ -1 ,
and by (3.17)

fx (χ) = R h w ≤ C χ -1 ,
which is the required conclusion. The proof for fy is entirely analogous.

Lemma 3.3 There exists a constant C such that

∆ h u -1 ≤ C ∇u , ∀u ∈ H 1 . (3.18) Proof: Let u ∈ H 1 . Then, for each w = 0 ∈ H 1 we have ( ∆ h u, w) w 1 = ( ∆ h u, P h w) w 1 ≤ ∇u ∇ P h w w 1 ≤ C ∇u ,
from which, taking the supremum over w ∈ H 1 we obtain (3.18). We used the stability of P h on H 1 , i.e.

the inequality

P h w 1 ≤ C w 1 , for w ∈ H 1 ,
which follows by the argument in [CT].

We now state and prove the main result of this section.

Theorem 3.1 For h sufficiently small, the semidiscrete problem (3.14) has a unique solution (η h , u h , v h ) in the interval [0, T ] of maximal existence of the solution (η, u, v) of the ibvp (R). Moreover, for some

constant C = C(η, u, v, T ) independent of h we have η -η h + u -u h + v -v h ≤ Ch r , and
η -η h 1 + u -u h 1 + v -v h 1 ≤ Ch r-1 , for each t ∈ [0, T ].
Proof: We suppose that for some constant M there holds that η ∞ ≤ M , u ∞ ≤ M and v ∞ ≤ M for 0 ≤ t ≤ T . Then, from (3.9b) for h sufficiently small

η 0 h ∞ ≤ η 0 h -η 0 ∞ + η 0 ∞ = R h η 0 -η 0 ∞ + η 0 ∞ ≤ Cγ(h) η 0 r,∞ + η 0 ∞ < 2M.
Similar estimates hold for u 0 h and v 0 h . The o.d.e. system (3.14) has a unique solution locally in t. By continuity, we may assume that there exists

t h ∈ (0, T ] such that u h ∞ ≤ 2M , η h ∞ ≤ 2M and v h ∞ ≤ 2M for all t ≤ t h .
We let now

ρ = η -R h η, θ = R h η -η h , τ = v -R h v, ζ = R h v -v h , σ = u -R h u, ξ = R h u -u h , so that θ ∈ S h , ζ, ξ ∈ S h and η -η h = ρ + θ, u -u h = σ + ξ, v -v h = τ + ζ. By (R) and (3.14) we have θ t = ĝx (σ + ξ) + ĝy (τ + ζ) + ĝx (uη -u h η h ) + ĝy (vη -v h η h ), (3.19)
ξ t = fx (θ + ρ) + 1 2 fx (u 2 ) -fx (u 2 h ) + fx (v 2 ) -fx (v 2 h ) + c fx (∆η -∆ h η h ), (3.20) ζ t = fy (τ + ζ) + 1 2 fy (u 2 ) -fy (u 2 h ) + fy (v 2 ) -fy (v 2 h ) + c fy (∆η -∆ h η h ).
(3.21)

The equation (3.20) may be written as

ξ t = fx (θ + ρ) + 1 2 fx (u(σ + ξ)) + fx ((σ + ξ)u h ) + fx (v(τ + ζ)) + fx ((τ + ζ)v h ) + c fx (∆η -∆ h η h ).
Taking L 2 -norms and using Lemma 3.1, (3.15), (3.12), (3.6a), (3.6b) and (3.18) we obtain, for 0 ≤ t ≤ t h ,

ξ t ≤ fx (θ + ρ) + 1 2 fx (u(σ + ξ)) + fx ((σ + ξ)u h ) + fx (v(τ + ζ)) + fx ((τ + ζ)v h ) + c fx (∆η -∆ h η h ) ≤ C( θ + ρ + u(σ + ξ) + (σ + ξ)u h + v(τ + ζ) + (τ + ζ)v h + fx (∆η -∆ h R h η) + fx ( ∆ h ( R h η -η h ) ) ≤ C[ θ + ρ + u L ∞ ( σ + ξ ) + v L ∞ ( τ + ζ ) + v h L ∞ ( τ + ζ ) + u h L ∞ ( σ + ξ ) + ∆η -∆ h R h η + ∆ h ( R h η -η h ) -1 ] ≤ C h r + θ + ξ + ζ + ∆η -P h ∆η + R h η -P h η + R h η -η h 1 ,
and so

ξ t ≤ C(h r + θ 1 + ξ + ζ ). (3.22)
Similarly, we have for 0

≤ t ≤ t h ζ t ≤ C(h r + θ 1 + ξ + ζ ). (3.23)
The equation (3.19) may be written as

θ t = ĝx (σ + ξ) + ĝy (τ + ζ) + ĝx (u(ρ + θ)) + ĝx ((σ + ξ)η h ) + ĝy (v(ρ + θ)) + ĝy ((τ + ζ)η h ).
Hence, taking H 1 -norms and using Lemma 3.1, and (3.6a,b) we have for 0

≤ t ≤ t h θ t 1 ≤ C( σ + ξ + τ + ζ + u(ρ + θ) + (σ + ξ)η h + v(ρ + θ) + (τ + ζ)η h ) ≤ C( σ + ξ + τ + ζ + u L ∞ ( ρ + θ ) + ( σ + ξ ) η h L ∞ + v L ∞ ( ρ + θ ) + ( τ + ζ ) η h L ∞ ) ≤ C( σ + ξ + τ + ζ + ρ + θ ),
and therefore that

θ t 1 ≤ C(h r + θ + ξ + ζ ), 0 ≤ t ≤ t h . (3.24) From (3.22)-(3.24) we see that for 0 ≤ t ≤ t h 1 2 d dt ( θ 2 1 + ξ 2 + ζ 2 ) ≤ Ch 2r + C( θ 2 1 + ξ 2 + ζ 2 ),
from which, using Gronwall's inequality, we have for 0

≤ t ≤ t h θ 1 + ξ + ζ ≤ Ch r . (3.25)
Hence, for t ≤ t h there holds

η -η h + u -u h + v -v h ≤ Ch r . (3.26)
Furthermore for t ≤ t h , we have, by (3.8), (3.9b) and (3.25)

η h -η L ∞ ≤ η h -R h η L ∞ + R h η -η L ∞ ≤ Ch -1 η h -R h η + Cγ(h) = Ch -1 θ 1 + Cγ(h) ≤ Ch r-1 Therefore, η h L ∞ ≤ η L ∞ + η h -η L ∞ ≤ Ch r-1 + M < 2M
for h sufficiently small. Similar estimates hold for u h , v h . These contradict the maximal property of t h and we conclude that we may take t h = T .

Hence (3.26) holds up to t = T , giving the desired optimal-rate L 2 -estimate. The O(h r-1 ) H 1 estimate follows easily by (3.6a,b) and (3.7).

It is worthwhile to note that temporal derivatives of arbitrary order of the semidiscrete solution (η h , u h , v h ) are bounded on [0, T ] by constants independent of h, as the following proposition shows.

Proposition 3.1 For h sufficiently small, let (η h , u h , v h ) be the solution of the semidiscrete problem

(3.14) for t ∈ [0, T ].
Then, for j = 0, 1, 2, 3, . . ., there exist constants C j independent of h, such that

max 0≤t≤T ∂ j t η h 1 + ∂ j t u h + ∂ j t v h ≤ C j .
(3.27)

Proof: From Theorem 3.1 we have, for some constant C 0 independent of h,

max 0≤t≤T ( η h 1 + u h + v h ) ≤ C 0 . (3.28)
Now, from (3.14), (ii) of Lemma 3.1, and (3.28), there follows for 0

≤ t ≤ T η ht 1 ≤ ĝx (u h ) 1 + ĝy (v h ) 1 + ĝx (η h u h ) 1 + ĝy (η h v h ) 1 ≤ C( u h + v h + η h u h + η h v h ) ≤ C(1 + η h ∞ ). (3.29)
In addition, from (3.14), (i) of Lemma 3.1, (3.15), and (3.28) we have for 0

≤ t ≤ T u ht ≤ fx (η h ) + 1 2 fx (u 2 h ) + 1 2 fx (v 2 h ) + |c| fx ( ∆ h η h ) ≤ C( η h + u 2 h + v 2 h + ∆ h η h -1 ) ≤ C(1 + u h ∞ u h + v h ∞ v h + η h 1 ) ≤ C(1 + u h ∞ + v h ∞ ).
(3.30)

A similar inequality holds for v ht . Now, from the closing argument of the proof of Theorem 3.1 we

may infer that max 0≤t≤T ( η h ∞ + u h ∞ + v h ∞ ) ≤ C, (3.31)
and, consequently, in view of (3.29) and (3.30), the validity of (3.27) for j = 1. Differentiating now the equations in (3.14) with respect to t and using again Lemma 3.1, (3.27) for j = 1, and (3.31) we see that (3.27) holds for j = 2 as well.

If we take the second temporal derivative of both sides of the first o.d.e. in (3.14), we see that in order to obtain a bound for ∂ 3 t η h 1 , we need, in addition to already established estimates, an h-independent bound for η h t ∞ . This is obtained as follows: For 0 ≤ t ≤ T we have, for θ = R h ηη h , from (3.8),

(3.9b), and (3.24), that

η h t ∞ ≤ θ t ∞ + R h η t ∞ ≤ Ch -1 θ t + Cγ(h) + C ≤ C.
We similarly get h-independent bounds for u ht ∞ and v ht ∞ that yield in turn similar bounds for ∂ 3 t u h and ∂ 3 t v h .

Hence (3.27) holds for j = 3 too. The case j = 4 follows immediately, as it does not need any L ∞ bounds on temporal derivatives of the semidiscrete approximations of order higher than one. To obtain (3.27) for j = 5 one needs, in addition to already established bounds, h-independent bounds for

∂ 2 t η h ∞ , ∂ 2 t u h ∞ ,
and

∂ 2 t v h ∞ on [0, T ].
These may be derived as follows: Differentiate with respect to t the expression for θ t after (3.23) and use the uniform bound on η ht ∞ and (3.22)-(3.25) to obtain that θ tt 1 ≤ Ch r .

The required bound for ∂ 2 t η h ∞ follows then from (3.8) and (3.9b). Similarly, differentiating e.g. the expression for ξ t after (3.21) and using the bounds on u ht ∞ , v ht ∞ and (3.22)-(3.25) we obtain that

ξ tt ≤ Ch r , from which ∂ 2 t u h ∞ ≤ C follows.
The case j = 6 requires no additional L ∞ bounds.

We continue by induction.

If j = 2k + 1, L ∞ -bounds for ∂ k t η h , ∂ k t u h , ∂ k t v h are found by differentiating
the expressions for θ t , ξ t and ζ t and using previously established bounds. The even case j = 2k + 2 requires no additional L ∞ bounds. (As a corrolary from the above proof it follows also that max 0≤t≤T

( ∂ j t η h ∞ + ∂ j t u h ∞ + ∂ j t v h ∞ ) ≤ C ′ j ,
holds for j = 0, 1, 2, . . ., where C ′ j are constants independent of h.)

From the result of this proposition we see that the system of o.d.e.'s (3.14) is not stiff. Therefore, one may use explicit time-stepping schemes to discretize (3.14) in the temporal variable without imposing stability mesh conditions on the time step ∆t in terms of h. Error estimates of optimal order in space and time in the case of explicit Runge-Kutta full discretizations may be established along the lines of the proof of Proposition 10 of [ADMII].

Remark 3.1 The H 1 error estimate in Theorem 3.1 may be strengthened as follows. For w ∈ H 1 define the discrete norm • 2,h as

w 2,h := ( w 2 1 + ∆ h w 2 ) 1/2 , ∀w ∈ H 1 .
Then, by letting w ∈ H 1 0 and considering the boundary-value problem f -b∆f = -w x in Ω with ∂f ∂n = 0 on ∂Ω, we easily see that ĝx (w) = R h f . A straightforward computation using (3.12) yields then that ĝx (w) 2 2,h ≤ C w x 2 , ĝy (w) 2 2,h ≤ C w y 2 . We may take now the • 2,h norm in (3.19) and obtain instead of (3.24). We conclude, along the lines of the proof of Theorem 3.1, that

θ t 2,h ≤ C(h r + θ + ξ + ζ ), 0 ≤ t ≤ t h ,
η -η h 2,h + u -u h 1 + v -v h 1 ≤ Ch r-1 .

Numerical experiments

In this section we present the results of numerical experiments that we performed in the case of the Bona-Smith systems using the modified Galerkin method as a base spatial discretization scheme. For the temporal discretization of the system of o.d.e.'s (3.14) we used an explicit, second-order Runge-Kutta method, the so-called "improved Euler" scheme. For the solution of the resulting linear systems at each time step we used the Jacobi-Conjugate Gradient method of ITPACK, taking the relative residuals equal to 10 -7 for terminating the iterations at each time step. In the computation of the discrete Laplacian ∆ h we used lumping of the mass matrix. In what follows we present numerical results that confirm the expected rates of convergence of the fully discrete scheme and three numerical experiments illustrating the use of the method in various surface flows of interest.

Experimental rates of convergence

In order to check the spatial convergence rates of the fully discrete modified Galerkin method we applied the scheme to the ibvp (R) for the Bona-Smith system with θ 2 = 9/11 taking as exact solution η(x, y, t) = cos(πx) cos(πy)e t , u(x, y, t) = x cos((πx)/2) sin(πy)e t , v(x, y, t) = y cos((πy)/2) sin(πx)e t .

defined on [0, 1] × [0, 1], which was covered by a uniform mesh consisting of isosceles right-angle triangles with perpendicular sides of length h = 2/N , where N is the number of the triangles. The solution was approximated in space by continuous, piecewise linear functions on this triangulation. We took ∆t = 0.01 and computed up to T = 1. In Tables 4.1 and 4.2 we show the resulting L 2 and H 1 errors and the corresponding experimental convergence rates that confirm the result of Theorem 3.1.

Table 4.1: L 2 errors and convergence rates for the modified Galerkin method for the Bona-Smith system with θ 2 = 9/11. Linear elements on triangular mesh and second-order RK time-stepping.

N η -η h rate(η) u -u h rate(u) v -v h rate(v) 512
1.3016E-2 -9.0816E-3 -9.0429E-2 -2048 3.2571E-3 1.9986 2.3603E-3 1.9440 2.3604E-2 1.9378 8192 8.1485E-4 1.9990 5.9882E-4 1.9788 5.9882E-4 1.9788 32768 2.0523E-4 1.9893 1.5417E-4 1.9576 1.5417E-4 1.9576 Table 4.2: H 1 errors and convergence rates for the modified Galerkin method for the Bona-Smith system with θ 2 = 9/11. Linear elements on triangular mesh and second-order RK time-stepping.

N η -η h 1 rate(η) u -u h 1 rate(u) v -v h 1 rate(v) 512 4.3739E-1 - 1.8597E-1 - 1.8488E-1 - 2048 2.
2030E-1 0.9894 8.7012E-2 1.0958 8.7008E-2 1.0874 8192 1.1038E-1 0.9970 4.2013E-2 1.0504 4.2013E-2 1.0503 32768 5.5221E-2 0.9992 2.0680E-2 1.0226 2.0680E-2 1.0226 cylinder, perpendicular to the impinging wave. In general, we did not observe great differences in the behaviour of the solutions of the two systems.

Evolution and reflection at the boundaries of a 'heap' of water

The sequence of plots of Figure 4.6 shows the temporal evolution of the free surface elevation of an initial Gaussian 'heap' of water with η 0 (x, y) = 2e -((x+40) 2 +(y+40) 2 )/5 , v 0 (x, y) = 0, as it collapses, forms radial outgoing riples and interacts with the reflective boundaries near a corner of the square [-80, 80] × [-80, 80]. Again, no major differences were observed between the solutions of the two systems (BBM-BBM and Bona-Smith, θ 2 = 9/11). The computation was effected with N = 84992 triangular elements and ∆t = 0.1. triangles and ∆t = 0.01 sec). For the impinging wave at t = 30 sec we estimated that A/h 0 ∼ = 0.014, λ/h 0 ∼ = 16.6, so that the Stokes number is approximately equal to 3.9, within the range of validity of the Boussinesq systems. The incoming wave hits the pier and the port boundary at x = 0, is reflected backwards and interacts with the remaining boundary.

  Figure 4.7 shows the corresponding one-dimensional plots (along x = -40) of the free surface elevation for both systems at four temporal instanses.

4. 4

 4 Line wave impinging on a 'port' structure In this experiment we integrated the BBM-BBM and the Bona-Smith (θ 2 = 9/11) systems in dimensional variables. Our domain represents part of a 'port' of depth h 0 = 50m consisting of the rectangle [-250, 250] × [0, 2000] minus the rectangular 'pier' [0, 100] × [0, 700]. (All distances in meters). Normal reflective boundary conditions are assumed to hold along the boundary of the pier and the intervals [-250, 0] and [100, 250] of the x-axis, while homogeneous Neumann boundary conditions have been imposed on η and v on the remaining parts of the boundary. An initial wave with η 0 (x, y) = Ae -y-1500 6000 , A = 1 m, and u 0 (x, y) = 0 m/sec, v 0 (x, y) = -1 2 (η 0 (x, y) -1 4 η 2 0 (x, y)) m/sec travels mainly towards the negative y direction shedding a dispersive tail behind. (Both systems were integrated with N = 77632

Figure 4 . 3 :Figure 4 . 7 :Figure 4

 43474 Figure 4.2: Experiment 4.2 Free surface elevation at four time instances. BBM-BBM vs. Bona-Smith (θ 2 = 9/11) systems.

  Let H s = H s (Ω), s ∈ R, denote the L 2 -based, real-valued, Sobolev classes on Ω and H 1 0 the subspace of H 1 whose elements have zero trace on ∂Ω. In the sequel, we shall denote by • and (•, •)

	.2)
	We now describe the class of boundary conditions on ∂Ω under which the problem (2.1)-(2.2) will
	be solved.

* This work was supported in part by a French-Greek scientific cooperation grant for the period 2006-08, funded jointly by EGIDE, France, and the General Secretariat of Research and Technology, Greece. D. Mitsotakis was also supported by Marie Curie Fellowship No. PIEF-GA-2008-219399 of the European Commission.

Solitary-wave-like pulse hitting a cylindrical obstacle

In our first experiment the domain that we consider is the rectangular channel [-15, 15] × [-30, 50] in which is placed a vertical impenetrable cylinder centered at (0, 10) with radius equal to 1.5. We pose normal reflective boundary conditions on the boundary of the cylinder and along the lines y = -30 and y = 50, and homogeneous Neumann boundary conditions for η and v on the lateral boundaries x = ±15.

As initial conditions we use the functions

where A = 0.2 and c s = 1 + A 2 , which represent a good approximation to a line solitary wave, [DMSI], of the BBM-BBM system. We integrate under these initial and boundary conditions the Bona-Smith system (1.1) with θ 2 = 9/11 and compare its evolving solution with that of the BBM-BBM system, i.e.

(1.1) with θ 2 = 2/3. The Bona-Smith system was discretized in space by the modified Galerkin method (amended in a straightforward way to handle the Neumann conditions on v on the lateral boundaries) with continuous, piecewise linear elements on a triangulation consisting of N = 290112 trangles; this mesh was fine enough for the 'convergence' of the numerical solution. The same truangulation was used to discretize in space the BBM-BBM system with the standard Galerkin method, [DMSI]. Both systems were discretized in time with the improved Euler method with ∆t = 0.01.

The line solitary-wave-like pulse is centered at y = -10 at t = 0. It propagates, mainly in the positive y direction, travelling with a speed c s of about 1.05, cf. Figure 4.2. The bulk of the solitary wave travels past the cylinder, while smaller amplitude scattered waves are produced by the interaction of the wave with the obstacle. Figure 4.3 shows, for both systems, contours of the elevation of the wave, superimposed on velocity vector plots, near the obstacle during the passage of the solitary wave. Figure 4.4 shows the free surface elevation for both systems as a function of y at x = 0 at several temporal instances. It is worthwhile to note that the maximum run-up at the extreme upstream point (x, y) = (0, 8.5) of the cylinder was equal to z = 0.335107 (achieved at t = 16.72) for the BBM-BBM system; the corresponding run-up for the Bona-Smith system was z = 0.319960 at t = 16.74. The analogous values for the run-up on the downstream point (0, 11.5) of the cylinder were z = 0.234310, t = 22.71 for the BBM-BBM system, and z = 0.230590, t = 22.68 for the Bona-Smith system.