Olivier Beaumont
email: olivier.beaumont@inria.fr

Loris Marchal
email: loris.marchal@ens-lyon.fr

Analysis of Dynamic Scheduling Strategies for Matrix Multiplication on Heterogeneous Platforms

Keywords: F.2.0 [Theory of Computation]: Analysis of Algorithms and Problem Complexity-General Algorithms, Performance Dynamic scheduling, data-aware algorithms, randomized algorithms, performance evaluation, matrix multiplication

The tremendous increase in the size and heterogeneity of supercomputers makes it very difficult to predict the performance of a scheduling algorithm. Therefore, dynamic solutions, where scheduling decisions are made at runtime have overpassed static allocation strategies. The simplicity and efficiency of dynamic schedulers such as Hadoop are a key of the success of the MapReduce framework. Dynamic schedulers such as StarPU, PaRSEC or StarSs are also developed for more constrained computations, e.g. task graphs coming from linear algebra. To make their decisions, these runtime systems make use of some static information, such as the distance of tasks to the critical path or the affinity between tasks and computing resources (CPU, GPU,. . . ) and of dynamic information, such as where input data are actually located. In this paper, we concentrate on two elementary linear algebra kernels, namely the outer product and the matrix multiplication. For each problem, we propose several dynamic strategies that can be used at runtime and we provide an analytic study of their theoretical performance. We prove that the theoretical analysis provides very good estimate of the amount of communications induced by a dynamic strategy and can be used in order to efficiently determine thresholds used in dynamic scheduler, thus enabling to choose among them for a given problem and architecture.

INTRODUCTION

Recently, there has been a very important change in both parallel platforms and parallel applications. On the one hand, computing platforms, either clouds or supercomputers involve more and more computing resources. This scale change poses many problems, mostly related to unpredictability and failures. Due to the size of the platforms, their complex network topologies, the use of heterogeneous resources, NUMA effects, the number of concurrent simultaneous computations and communications, it is impossible to predict exactly the time that a specific task will take. Unpredictability makes it impossible to statically allocate the tasks of a DAG onto the processing resources and dynamic scheduling and allocation strategies are needed. As a consequence, in recent years, there has been a large amount of practical work to develop efficient runtime schedulers. The main characteristics of these schedulers is that they make their decisions at runtime, based on the expected duration of the tasks on the different kind of processing units (CPUs, GPUs,...) and on the expected availability time of the task input data, given their actual locations. Thanks to these information, the scheduler allocates the task to the resource that will finish its processing as soon as possible. Moreover, all these runtime systems also make use of some static information that can be computed from the task graph itself, in order to decide the priority between several ready tasks. This information mostly deals with the estimated critical path as proposed in HEFT [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF] for instance.

On the other hand, there has been a dramatic simplification of the application model in many cases, as asserted by the success of the MapReduce framework [START_REF] Dean | MapReduce: Simplified data processing on large clusters[END_REF] which has been popularized by Google. It allows users without particular knowledge in parallel algorithms to harness the power of large parallel machines. In MapReduce, a large computation is broken into small tasks that run in parallel on multiple machines, and scales easily to very large clusters of inexpensive commodity computers. MapReduce is a very successful example of dynamic schedulers, as one of its crucial feature is its inherent capability of handling hardware failures and processing capabilities heterogeneity, thus hiding this complexity to the programmer, by relying on on-demand allocations and the on-line detection of nodes that perform poorly (in order to re-assign tasks that slow down the process). As we explained in a previous work [START_REF] Beaumont | Non linear divisible loads: There is no free lunch[END_REF], MapReduce, although tailored for linear complexity operations (such as text parsing), is now widely used for non linear complexity tasks. In this case, it induces a large replication of the data.

For example, when MapReduce is used to compute the outer product of two vectors a and b, the most common technique is to emit all possible pairs of (ai, bj ), so that many processors can be used to compute the elementary products. This induces a large replication factor, since MapReduce is not aware of the 2-dimensional nature of the data.

Our goal in this paper is to show how simple data-aware dynamic schedulers can be proven efficient in a specific context. We concentrate here on two elementary kernels, namely the outer product and the matrix multiplication. These kernels do not induce dependencies among their tasks, but because of their massive input data reuse results, a straightforward MapReduce implementations of these kernels would involve a large replication overhead. Indeed, in both cases [START_REF] Beaumont | Non linear divisible loads: There is no free lunch[END_REF], input vectors or input matrices need to be replicated when the kernel is processed by a large-scale parallel platform, and basic dynamic strategies that allocate tasks at random to processors fail to achieve reasonable communication volumes with respect to known lower bounds.

In the present paper, we first present and study a very simple yet efficient dynamic scheduler for the outer product, that generates a communication volume close to the lower bound. Our main contribution is to analyze the communication volume generated by the dynamic scheduler as a continuous process that can be modeled by an Ordinary Differential Equation (ODE). We prove that the analytic communication volume of the solution of the ODE is close to the actual communication volume as measured using simulations. Moreover, we prove that this analysis of the solution of the ODE can be used in order to optimize a dynamic randomized allocation strategy, for instance, by switching between two strategies when the number of remaining tasks is smaller than a given threshold, that is determined by the theoretical analysis. This simple example attests the practical interest of the theoretical analysis of dynamic schedulers, since it shows that the analytic solution can be used in order to incorporate static knowledge into the scheduler. After presenting our method on the outer product (Section 3), we move to a more common kernel, the matrix multiplication and show how the previous analysis can be extended in Section 4.

RELATED WORK

We briefly review previous works related to our study, which deals both with actual runtime schedulers and with their theoretical studies.

Runtime dynamic schedulers

As mentioned in the introduction, several runtime systems have been recently proposed to schedule applications on parallel systems. Among other successful projects, we may cite StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF], from INRIA Bordeaux (France), DAGuE and PaRSEC [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for high performance computing[END_REF][START_REF] Bosilca | PaRSEC: A programming paradigm exploiting heterogeneity for enhancing scalability[END_REF] from ICL, Univ. of Tennessee Knoxville (USA) StarSs [START_REF] Planas | Hierarchical task-based programming with StarSs[END_REF] from Barcelona Supercomputing Center (Spain) or KAAPI [START_REF] Gautier | Kaapi: A thread scheduling runtime system for data flow computations on cluster of multi-processors[END_REF] from INRIA Grenoble (France). Most of these tools enable, to a certain extent, to schedule an application described as a task graph (usually available in the beginning of the computation, but sometimes generated and discovered during the execution itself), onto a parallel platforms. Most of these tools allow to harness complex platforms, such as multicores and hybrid platforms, including GPUs or other accelerators. These runtime systems usually keep track of the occupation of each computing devices and allocate new tasks on the processing unit that is expected to minimize its completion time. Our goal in this paper in to provide an analysis of such dynamic schedulers for simple operations, that do not involve tasks dependencies but massive data reuse.

Theoretical studies of dynamic systems

Many studies have proposed to use queuing theory [START_REF] Gross | Fundamentals of Queueing Theory[END_REF] to study the behavior of simple parallel systems and their dynamic evolution. Among many others, Berten et al. [START_REF] Berten | Brokering strategies in computational grids using stochastic prediction models[END_REF] propose to use such stochastic models in order to model computing Grids, and Mitzenmacher [START_REF] Mitzenmacher | How useful is old information?[END_REF] studies how notto-date information can lead to bad scheduling decisions in a simple parallel system.

Recently, mean field techniques [START_REF] Gast | Mean field for Markov decision processes: from discrete to continuous optimization[END_REF][START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF] have been proposed for analyzing such dynamic processes. They give a formal framework to derive a system of ordinary differential equations that is the limit of a Markovian system when the number of objects goes to infinity. Such techniques have been used for the first time in [START_REF] Mitzenmacher | Analyses of load stealing models based on differential equations[END_REF] where the author derives differential equations for a system of homogeneous processors who steal a single job when idle.

RANDOMIZED DYNAMIC STRATEGIES FOR THE OUTER-PRODUCT

We present here the analysis of a dynamic scheduler for a simple problem from linear algebra, namely the outerproduct of two vectors.

Problem definition

We consider the problem of computing the outer-product ab t of two large vectors a and b of size N , i.e. to compute all values ai × bj, ∀1 ≤ i, j ≤ N . The computing domain can therefore be seen as a matrix. For granularity reasons, we will consider that a and b are in fact split into N/l blocks of size l and that a basic operation consists in computing the outer product of two (small) vectors of size l.

As stated above, we target heterogeneous platforms consisting of p processors P1, . . . , Pp, where the speed of processor Pi, i.e. the number of outer products of size l vectors that P k can do in one time unit, is given by s k . We will also denote by rs k the relative speed of rs k = s k i s i . Note that the randomized strategies that we propose are agnostic to processor speeds, but they are demand driven, so that a processor with a twice larger speed will request work twice faster.

In the following, we assume that a master processor coordinates the work distribution: it is aware of which a and b blocks are replicated on the computing nodes and decides which new blocks are sent, as well as which tasks are allocated to the nodes. After completion of their allocated tasks, computing nodes simply report to the master processor, requesting for new tasks.

We will assume throughout the analysis that it is possible to overlap computations and communications. This can be achieved with dynamic strategies by uploading a few blocks in advance at the beginning of the computations and then to request work as soon as the number of blocks to be processed becomes smaller than a given threshold. Determining this threshold would require to introduce a communication model and a topology, what is out of the scope of this paper, and we will assume that the threshold is known. In practice, the number of tasks required to ensure a good overlap has been observed to be small in [START_REF] Kreaseck | Autonomous protocols for bandwidth-centric scheduling of independent-task applications[END_REF][START_REF] Parashar | Autonomic computing: concepts, infrastructure, and applications[END_REF] even though a rigorous algorithm to estimate it is still missing.

As we observed [START_REF] Beaumont | Non linear divisible loads: There is no free lunch[END_REF], performing a non linear complexity task such as a Divisible Load or a MapReduce operation requires to replicate initial data. Our objective is to minimize the overall amount of communications, i.e. the total amount of data (the number of blocks of a and b) sent by the master node initially holding the data, or equivalently by the set of devices holding the data since we are interested in the overall volume only, under the constraint that a perfect loadbalancing should be achieved among resources allocated to the outer product computation. Indeed, due to data dependencies, if we were to minimize communications without this load-balancing constraint, the optimal (but very inefficient) solution would consist in making use of a single computing resource so that each data block would be sent exactly once.

Design of randomized dynamic strategies

As mentioned above, vectors a and b are split into N/l data blocks. In the following, we denote by ai the ith block of a (rather than the ith element of a) since we always consider elements by blocks. As soon as a processor has received two data blocks ai and bj, it can compute the block Mi,j = (ab t )i,j = aib t j . This elementary task is denoted by Ti,j . All data blocks are initially available at the master node only.

One of the simplest strategy to allocate computational tasks to processors is to distribute tasks at random: whenever a processor is ready, a task Ti,j is chosen uniformly at random among all available tasks and is allocated to the processor. The data corresponding to this task that is not yet on the processor, that is one or two of the ai and bj blocks are sent by the master. We denote this strategy by Random-Outer. Another simple option is to allocate tasks in lexicographical order of indices (i, j) rather than randomly. This strategy will be denoted as SortedOuter.

Both previous algorithms are expected to induce a large amount of communications because of data replication. Indeed, in these algorithms, there is no reason why the data sent for the processing of tasks on a given processor P k may be re-used for upcoming tasks. This is why dynamic dataaware strategies have been introduced. In the runtime systems cited above, such as StarPU, the scheduler is aware of the locality of the data and uses this information when allocating tasks to processors: it is much more beneficial, when allocating a new task on P k , to take advantage of the a and b data already present on the processor, and to compute for example all possible products aib t j ′ before sending new blocks of data. We propose such a strategy, denoted DynamicOuter, in Algorithm 1: when a processor P k receives a new pair of blocks (ai, bj ), all possible products aib t j ′ and a i ′ b t j are also allocated to P k , for all data blocks a i ′ and b j ′ that have already been transmitted to P k in previous steps.

Note that the DynamicOuter scheduler is not computationally expensive: it is sufficient to maintain a set of unknown a and b data (of size O(N/l)) for each processor, and to randomly pick an element of this set when allocating new blocks to a processor P k .

We have compared the performance of previous schedulers through simulations on Figure 1. Processor speeds are chosen uniformly in the interval [START_REF] Gautier | Kaapi: A thread scheduling runtime system for data flow computations on cluster of multi-processors[END_REF]100], which means a large Algorithm 1: DynamicOuter strategy.

while there are unprocessed tasks do Wait for a processor P k to finish its tasks I ← {i such that P k owns ai} J ← {j such that P k owns bj} Choose i / ∈ I and j / ∈ J uniformly at random Send ai and bj to P k Allocate all tasks of {Ti,j } ∪ {T i,j ′ , j ′ ∈ J} ∪ {T i ′ ,j , i ′ ∈ I} that are not yet processed to P k and mark them processed

LB = 2N k √ rs k = 2N k s k i si
, where s k is the speed of processor P k and rs k its relative speed.

Indeed, in a very optimistic setting, each processor is dedicated to computing a "square" area of M = ab t , whose area is proportional to its relative speed, so that all processors finish their work at the same instant. In this situation, the amount of communications for P k is proportional to the half perimeter of this square of area N 2 rs k . Note that this lower bound is not expected to be achievable (consider for instance the case of 2 heterogeneous processors). The best known static algorithm (based on a complete knowledge of all relative speeds) has an approximation ratio of 7/4 [START_REF] Beaumont | Partitioning a square into rectangles: NP-completeness and approximation algorithms[END_REF]. This algorithm computes an allocation scheme based on the computing speeds of the processors. As outlined in the introduction, such an allocation mechanism is not practical in our context, since our aim is to rely on more dynamic runtime strategies, but can be used as a comparison basis.

As expected, we notice on Figure 1 that data-aware strategies induce significantly less communication than purely random strategies.

Our DynamicOuter allocation scheme suffers some limitation: when the number of remaining blocks to compute is small, the proposed strategy is inefficient as it may send a large number of a and b blocks to a processor P k before it is able to process one of the last few available tasks. Thus, we propose an improved version DynamicOuter2Phases in Algorithm 2: when the number of remaining tasks becomes smaller than a given threshold, we switch to the basic randomized strategy: any available task Ti,j is allocated to a requesting processor, without taking data locality into account. The corresponding data ai and bj are then sent to P k if needed.

Algorithm 2: DynamicOuter2Phases strategy. while the number of processors is larger than the threshold do

Wait for a processor P k to finish its tasks I ← {i such that P k owns ai} J ← {j such that P k owns bj} Choose i / ∈ I and j / ∈ J uniformly at random Send ai and bj to P k Allocate all tasks of {Ti,j } ∪ {T i,j ′ , j ′ ∈ J} ∪ {T i ′ ,j , i ′ ∈ I} that are not yet processed to P k and mark them processed while there are unprocessed tasks do Wait for a processor P k to finish its tasks Choose randomly an unprocessed task Ti,j if P k does not hold ai then send ai to P k if P k does not hold bj then send bj to P k Allocate Ti,j to P k As illustrated on Figure 2, for a well chosen number of tasks processed in the second phase, this new strategy allows to reduce further the amount of communications. However, this requires to accurately set the threshold, depending on the size of the matrix and the relative speed of the processors. If too many tasks are processed in the second phase, the performance is close to the one of RandomOuter. On the contrary, if too few tasks are processed in the second phase, the behavior becomes close to DynamicOuter. The optimal threshold corresponds here to a few percent of tasks being processed in the second phase. In the following, we present an analysis of the DynamicOuter2Phases strategy that both allows to predict its performance and to optimally set the threshold, so as to minimize the amount of communications.

Theoretical analysis of dynamic randomized strategies

In this section, our aim is to provide an analytical model for Algorithm DynamicOuter2Phases. Analyzing such a strategy is crucial in order to assess the efficiency of runtime dynamic strategies and in order to tune the parameters of dynamic strategies or to choose among different strategies depending on input parameters.

In what follows, we assume that N , the size of vectors a and b, is large and we consider a continuous dynamic process whose behavior is expected to be close to the one of DynamicOuter2Phases. In what follows, we concentrate on processor P k whose speed is s k . At each step, Dynamic-Outer2Phases chooses to send one data block of a and one data block of b, so that P k knows the same number of data blocks of a and b. As previously, we denote by M = ab t the result of the outer product and by Ti,j the tasks that corresponds to the product of data blocks ai and bj We denote by x = y/N the ratio of elements of a and b that are known by P k at a given time step of the process and by t k (x) the corresponding time step. We concentrate on a basic step of DynamicOuter2Phases during which the fraction of data blocks of both a and b known by P k goes from x to x + δx. In fact, since DynamicOuter2Phases is a discrete process and the ratio known by P k goes from x = y/N to x + l/N = y/N + l/N . Under the assumption that N is large, we assume that we can approximate the randomized discrete process by the continuous process described by the corresponding Ordinary Differential Equation on expected values. The proof of convergence is out of the scope of this paper but we will show that this assumption provides very good results through simulations in Section 3.4.

Let us remark that during the execution of Dynamic-Outer2Phases, tasks Ti,j are greedily computed as soon as a processor knows the corresponding data blocks of ai and bj. Therefore, at time t k (x), all tasks Ti,j such that P k knows data blocks ai and bj have been processed and there are x 2 N 2 /l 2 such tasks. Note also that those tasks may have been processed either by P k or by another processor Pj since processors compete to process tasks. Indeed, since data blocks of a and b are possibly replicated on several processors, then both P k and P ℓ may know at some point both ai and bj . In practice, the processor which computes Ti,j is the one that learns both ai and bj first.

Figure 3 depicts the computational domain during the first phase of DynamicOuter2Phases from the point of view of a given processor P k (rows and columns have been reordered for the sake of clarity). The top-left square (in blue) corresponds to value of a and b that are known by P k , and all corresponding tasks have already been processed (either by P k or by another processor). The remaining "L"-shaped area (in grey) corresponds to tasks Ti,j such that P k does not hold either the corresponding value of a, or the corresponding value of b, or both. When receiving a new value of a and b (in red), P k is able to process all the tasks (in red) from the two corresponding row and column. Some elements The top-left blue rectangle represents the data owned by the processor at time t k (x) (a permutation of the rows and columns has been applied to have it in the upper left corner). The new elements δx are depicted in red, as well as the corresponding available tasks. Note that some tasks (in black) corresponding to the combination of δx with the known elements have already been processed by other processors.

from this row and this column may be already processed (in black).

In what follows, we denote by g k (x) the fraction of tasks Ti,j in the previously described "L"-shaped area that have not been computed yet. We also assume that the distribution of unprocessed tasks in this area is uniform, and we claim that this assumption is valid for a reasonably large number of processors. Our simulations below show that this leads to a very good accuracy.

Based on this remark, we are able to prove the following Lemma

Lemma 1. g k (x) = (1 -x 2 ) α k , where α k = i =k s i s k .
Proof. Let us consider the tasks that have been computed by all processors between t k (x) and t k (x + δx). As depicted on Figure 3, these tasks can be split into two sets.

• The first set of tasks consists in those that can be newly processed by P k between t k (x) and t k (x + δx). P k has the possibility to combine the δxN new elements of a with the xN already known elements of b (and to combine the δxN new elements of b with the xN already known elements of a). There is therefore a total of 2 x δx N 2 such tasks (at first order). Among those, by definition of g, the expected number of tasks that have not already been processed by other processors is given by 2 x δx g(x) N 2 . Therefore, the expected duration of this step is given by t

k (x+δx)-t k (x) = 2 x δx g(x) N 2 s k .
• The second set of tasks consists in those computed by other processors Pi, i = k. Our assumption states that we are able to overlap communications by computations (by uploading data blocks slightly in advance), so that processors Pi, i = k will keep processing tasks between t k (x) and t k (x + δx) and will process on expectation 2 x δx g(x) N 2 i =k s i s k tasks.

Therefore, we are able to estimate how many tasks will be processed between t k (x) and t k (x+δx) and therefore to compute the evolution (on expectation) of g k . More specifically, we have

g k (x + δx) 1 -(x + δx) 2 N 2 = g k (x) (1-x 2 ) N 2 -2 x δx g(x) N 2 -2 x δx g(x) N 2 i =k si s k ,
which gives at first order

g k (x + δx) -g k (x) = g k (x) δx -2 x α k 1 -x 2 ,
where

α k = i =k s i s k
. Therefore, the evolution of g k with x is given by the following ordinary differential equation

g ′ k (x) g k (x) = -2 x α k 1 -x 2
where both left and right terms are of the form

f ′ /f , what leads to ln(g k (x)) = α k ln(1 -x 2 ) + K
and finally to

g k (x) = exp(K)(1 -x 2 ) α k ,
where exp(K) = 1 since g k (0) = 1. This achieves the proof of Lemma 1.

Remember that t k (x) denotes the time necessary for P k to know x elements of a and b. Then,

Lemma 2. t k (x) i si = N 2 (1 -(1 -x 2 ) α k +1 ).
Proof. We have seen that some of the tasks that could have been processed by P k (tasks Ti,j such that P k knows both ai and bj) have indeed been processed by other processors. In order to prove the lemma, let us denote by h k (x) the number of such tasks at time t k (x). Then

h k (x + δx) = h k (x) + 2 x δx (1 -g k (x))N 2
by definition of g k so that, using Lemma 1,

h ′ k (x) = N 2 (2x -2x(1 -x 2 ) α k ) and h k (x) = N 2 (x 2 + (1 -x 2 ) α k +1 α k + 1 + K)
and since h k (0) = 0,

h k (x) = N 2 (x 2 + (1 -x 2 ) α k +1 α k + 1 - 1 α k + 1
).

Moreover, at time t k (x), all the tasks that could have been processed by P k have • either been processed by P k and there are exactly t k (x)s k such tasks since P k has been processing all the time in this area,

• or processed by other processors and there are exactly h k (x) such tasks by definition of h k .

Therefore,

x 2 N 2 = h k (x) + t k (x)s k and finally t k (x) i si = N 2 (1 -(1 -x 2 ) α k +1 ),
which achieves the proof of Lemma 2.

Above equations well describe the dynamics of Dynamic-Outer2Phases as long as it is possible to find blocks of a and b that enable to compute enough unprocessed tasks. On the other hand, at the end, it is better to switch to another algorithm, where unprocessed tasks Ti,j are picked up randomly, which possibly requires to send two blocks ai and bj. In order to decide when to switch from one strategy to the other, we introduce an additional parameter β.

As presented above, a lower bound on the communication volume received by P k (if perfect load balancing is achieved) is given by LB = 2N k √ rs k . We will switch from the DynamicOuter to the RandomOuter strategy when the fraction of tasks x 2 k N 2 for which P k owns the input data is approximately β times what it would have computed optimally, that is, when x 2 k is close to β s k i s i = βrs k , for a value of β that is to be determined. For the sake of the analysis, it is important that we globally define the instant at which we switch to the random strategy, and that it does not depend on the processor P k . In order to achieve this, we look for x 2 k as x 2 k = (βrs k -αrs 2 k ) and we search α such that t k (x k ) does not depend on k at first order in 1/rs k , where rs k is of order 1/p and p is the number of processors.

Lemma 3. If α = β 2 /2, then t k (x k ) i si = N 2 (1 -e -β (1 + o(rs k ))). Proof. Since t k (x k ) i si = N 2 (1 -(1 -x 2 k ) α k +1 , then t k (x k ) = N 2 i si (1 -e 1 rs k ln(1 -βrs k + αrs 2 k )) = N 2 i si (1 -e 1 rs k (-βrs k + αrs 2 k -(βrs k ) 2 /2)) (at first order) = N 2 i si (1 -e -β (1 + o(rs k ))).
which achieves the proof of Lemma 3.

One remarkable characteristics of the above result is that it does not depend (at least up to order 2) on k anymore. Otherwise stated, at time T = N 2 i s i (1 -e β ), each processor P k has received (βrs k -

β 2 /2rs 2 k )N 2 = √ βrs k (1 -
βrs k /4)N data, to be compared with the lower bound on communications for processor P k : √ rs k N .

Using both these results, it is possible to derive the ratio between the overall amount of communication induced by the first phase with respect to the lower bound as a function of β.

Lemma 4. Let us denote by VPhase1 the volume of the communications induced by Phase 1 and by LB = 2N k

√ rs k the lower bound for the communications induced by the whole outer product, then

VPhase1 LB ≤ β + β 3/2 i rs 3/2 k

4LB

(at first order).

Proof. The proof is obtained by replacing VPhase1 by k √ βrs k (1 -βrs k /4)N .

Lemma 4 provides the evaluation of the expected communication volume induced by the first phase of Dynamic-Outer2Phases with respect to the lower bound. In the following, we will establish a similar result for the second phase in Lemma 5.

Lemma 5. Let us denote by VPhase2 the volume of the communications induced by Phase 1 and by LB = 2N k

√ rs k the lower bound for the communications induced by the whole outer product, then

VPhase2 LB ≤ e -β N 1 - √ β k rs 3/2 k k rs 1/2 k
(at first order).

Proof. During Phase 2, when a processor P k requests some work, a random task is sent among those that have not been processed yet. This task Ti,j induces either the communication of one data block (if either ai or bj is already know at P k ) or 2 data blocks (but not 0 by construction).

More precisely, since tasks are sent at random and since P k knows a fraction x k = √ βrs k (1-βrs k /4) of the elements of a and b at the end of Phase 1,

• a task induces the communication of one block with probability 2x k 1+x k ,

• a task induces the communication of two blocks with probability 1-x k 1+x k . so that the expected number of communications per task for

P k is 2x k 1 + x k × 1 + 1 -x k 1 + x k × 2 = 2 1 + x k .
Moreover, since Phase 2 starts at the same instant on all processors and since processors are continuously processing tasks, P k processes a fraction rs k of the e -β N 2 remaining tasks. The overall communication cost induced by Phase 2 is therefore given (on expectation and at first order) by

VPhase2 = e -β N 2 1 -β k rs 3/2 k
, which achieves the proof of Lemma 5.

Theorem 6. The ratio of the overall volume of communications to the lower bound if we switch between both phases when e -β N 2 tasks remain to be processed is given by

β + β 3/2 k rs 3/2 k 4 k rs 1/2 k + e -β N 2 1 - √ β k rs 3/2 k k rs 1/2 k .
Theorem 6 is a direct consequence of Lemma 4 and Lemma 5. Therefore, in order to minimize the overall amount of communications, we numerically determine the value of β that minimizes the above expression and then switch between Phases 1 and 2 when e -β N 2 tasks remain to be processed.

Assessing the validity of the analysis through simulations

We have performed simulations to study the accuracy of the previous theoretical analysis, that is a priori valid only for large values of p and N/l, and to show how it is helpful to compute the threshold for DynamicOuter2Phases. The simulations have been done using an ad-hoc event based simulation tool, where processors request new tasks as soon as they are available, and tasks are allocated based on the given runtime dynamic strategy. Again, processor speeds are chosen uniformly in the interval [START_REF] Gautier | Kaapi: A thread scheduling runtime system for data flow computations on cluster of multi-processors[END_REF]100]. This degree of heterogeneity may seem excessive but we show in Section 3.5 that using a different heterogeneity model does not significantly impact the results. The communication amount of each strategy is normalized by the lower bound computed in Section 3.3. Figure 4 presents the results for vectors of 100 blocks and Figure 5 does the same for vectors of 1000 blocks. In both figures, the analysis is extremely close to the performance of DynamicOuter2Phases (which makes them indistinguishable on the figures) and proves that our analysis succeed to accurately model our dynamic strategy, even for relatively small values of p and N/l. Moreover, we can see in Figure 5 that it is even more crucial to use a dataaware dynamic scheduler when N is large, as the ratio between the communication amount of simple random strategies (RandomOuter and SortedOuter) and dynamic dataaware schedulers (such as DynamicOuter2Phases) can be very large.

Our second objective is to show that the theoretical analysis that we propose can be used in order to accurately compute the threshold of DynamicOuter2Phases, i.e., that the β parameter computed earlier is close to the best one. To do this, we compare the communication amount of Dynamic-Outer2Phases for various values of the β parameter. Figure 6 shows the results for 20 processors and N/l = 100. This is done for a single and arbitrary distribution of computing speeds, as it would make no sense to compute average values for different distributions since they would lead to different optimal values of β. This explains the irregular performance graph for DynamicOuter2Phases. This figure shows that in the domain of interest, i.e. for 3 ≤ β ≤ 6, the analysis correctly fits to the simulations, and that the value of β that minimizes the analysis (here β = 4.17) lies the interval of β values that minimize the communication amount of DynamicOuter2Phases. To compare to Figure 2, this corresponds to 98.5% of the tasks being processed in the first phase.

Impact of the heterogeneity

The speed distribution used in the previous experiments (speeds taken in the interval [START_REF] Gautier | Kaapi: A thread scheduling runtime system for data flow computations on cluster of multi-processors[END_REF]100]) may seem too heterogeneous to reasonably model actual computing platforms, where heterogeneity comes either from the use of a few classes of different processors (new and old machines, processor equipped with accelerators or not, etc.) or from the fact that machines are not dedicated, which implies stochastically variable processor speed. It is natural to ask whether the speed distribution impacts the ranking of the previous heuristics, or the accuracy of our analysis. Figure 7 presents the behavior of all previous heuristics for a varying range of heterogeneity. A heterogeneity of 0 means perfectly homogeneous computing speeds, while a heterogeneity of 100 means that the ratio between the smallest and the largest speeds is large. In this figure and the following one, error bars represents the standard deviations with 50 tries. We notice that the heterogeneity degree has very little impact on the relative amounts of communication of the studied heuristics. In Figure 8, we study the same heuristics using different scenarios:

• Scenarios unif.1 and unif.2 corresponds to the previous setting, with speeds taken uniformly at random in intervals [80, 120] (unif.1) and [50, 150] (unif.2).

• Scenarios set.3 and set.5 corresponds to the case when there are a few classes of processors with different speed. The speeds are then taken uniformly from the set of possible speeds: (80, 100, 150) for set.3 or (40, 80, 100, 150, 200) for set.5.

• Scenarios dyn.5 and dyn.20 corresponds to very simple dynamic settings. Each computing speed is first taken uniformly at random in interval [80,120]. Then, after computing a task, a processor sees its computing speed randomly changed by up to 5% (dyn.5) or 20% (dyn.20).

This figure shows that neither the speed distribution nor the dynamic evolution of the speeds notably affect the performance of the heuristics.

Runtime estimation of β

In order to estimate the β parameter in the Dynamic-Outer2Phases strategy, it seems necessary to know the processing speed, as β depends on k s k / i si. However, we have noticed a very small deviation of β with the speeds. For example, in Figure 6, the value of β computed when assuming homogeneous speeds (4.1705) is very close to the one computed for heterogeneous speeds (4.1679).

For a large range of N/l and p values (namely, p in [10, 1000] and N/l ∈ [max(10, √ p), 1000]), for processor speeds in [START_REF] Gautier | Kaapi: A thread scheduling runtime system for data flow computations on cluster of multi-processors[END_REF]100],

the optimal value for β goes from 1 to 6.2. However, for fixed values of N/l and p, the deviations among the β values obtained for different speed distributions is at most 0.045 (with 100 tries). Our idea is to approximate β with β hom computed using a homogeneous platform with the same number of processors and with the same matrix size. The relative difference between β hom and the average β of the previous set is always smaller than 5%. Moreover, the error on the communication volume predicted by the analysis when using homogeneous speeds instead of the actual ones is at most 0.1%. These figures are derived with the most heterogeneous speed distribution (speeds in [START_REF] Gautier | Kaapi: A thread scheduling runtime system for data flow computations on cluster of multi-processors[END_REF]100]) and thus hold for the other distributions of Section 3.5 as well. This proves that even if our previous analysis ends up with a formula for β that depends on the computing speeds, in practice, only the knowledge of the matrix size and of the number of processors are actually needed to define the threshold β. Our dynamic scheduler DynamicOuter2Phases is thus totally agnostic to processor speeds.

MATRIX MULTIPLICATION

We adapt here the previous dynamic algorithm and its theoretical analysis to a more complex problem: the multiplication of two matrices.

Notations and dynamic strategies

We first adapt the notations to the problem of computing the product of two matrices C = AB. As in the previous section, we consider that all transfers and computations are performed using blocks of size l ×l, so that all three matrices are composed of N 2 /l 2 blocks and Ai,j denotes the block of A on the ith row and the jth column. The basic computation step is a task T i,j,k , which corresponds to the update Ci,j ← Ci,j + A i,k B k,j . To perform such a task, a processor has to receive the input data from A and B (of size 2l 2 ), and to send the result (of size (N/l) 2 ) back to the master at the end of the computation. Thus, it results in a total amount of communication of 3 (N/l) 2 . As previously, in order to minimize the amount of communications, our goal is to take advantage of the blocks of A, B and C that have already been sent to a processor Pu when allocating a new task to Pu. Note that at the end of the computation, all Ci,js are sent back to the master that computes in turn the final results by adding the different contributions. This computational load is much smaller than computing the products T i,j,k and we will neglect it in what follows.

As we assume that processors work during the whole process, the load imbalance, i.e. the difference between the amount of work processed by Pi and what it should have processed given its speed is at most one block. Thus, a maximal block size l can easily be derived from a maximal load imbalance. The value of l must also be large enough to overlap communications of size 3l 2 with computations of size l 3 . As usual, the block size should also be large enough to benefit from BLAS effect and small enough so as to fit into caches. We assume that the optimal block size l is computed by the runtime environment.

The simple strategies RandomOuter and SortedOuter translate very easily for matrix multiplication into the strategies RandomMatrix and SortedMatrix. We adapt the DynamicOuter strategy into DynamicMatrix as follows.

We ensure that at each step, for each processor Pu there exist sets of indices I, J and K such that Pu owns all values A i,k , B k,j , Ci,j for i ∈ I, j ∈ J and k ∈ K, so that it is able to compute all corresponding tasks T i,j,k . When a processor becomes idle, instead of sending a single block of A, B and C, we choose a tuple (i, j, k) of new indices (with i / ∈ I , j / ∈ J and k / ∈ K) and send to Pu all the data needed to extend the sets I, J, K with (i, j, k). This corresponds to sending 3 × (2|I| + 1) data blocks to Pu (note that |I| = |J| = |K|). In fact, blocks of C are not send by the master to the processor, but on the contrary will be sent back to the master at the end of the computation; however, this does not change the analysis since we are only interested in the overall volume of communications. Then, processor Pu is allocated all the unprocessed tasks that can be done with the new data. Algorithm 3 details this strategy.

As in the case of the outer product, when the number of remaining blocks to be processed becomes small,Random-Matrix strategy outperforms the DynamicMatrix strategy. Therefore, we introduce the intermediate Dynamic-Matrix2Phases strategy that consists into two phases. During Phase 1, the DynamicMatrix strategy is used. Then, when the number of remaining tasks becomes smaller than e -β N 3 for a value of β that is to be determined, we switch to Phase 2 and use strategy RandomMatrix. As in the case of the outer product, the theoretical analysis proposed in the next section will help us to determine the optimal value of β, i.e. the instant when to switch between phases in order to minimize the overall communication volume in the DynamicMatrix2Phases strategy.

Theoretical analysis of dynamic randomized strategies

In this section, our aim is to provide an analytical model for Algorithm DynamicMatrix2Phases similarly to what Algorithm 3: DynamicMatrix strategy.

while there are unprocessed tasks do Wait for a processor Pu to finish its task I ← {i such that Pu owns A i,k for some k} J ← {i such that Pu owns B k,j for some k} K ← {i such that Pu owns A i,k for some i} Choose i / ∈ I , j / ∈ J and k / ∈ K uniformly at random Send the following data blocks to Pu:

• A i,k ′ for k ′ ∈ K ∪ {k} and A i ′ ,k for i ′ ∈ I ∪ {i} • B k,j ′ for j ′ ∈ J ∪ {j} and B k ′ ,j for k ′ ∈ K ∪ {k} • C i,j ′ for j ′ ∈ J ∪ {j} and C i ′ ,j for i ′ ∈ I ∪ {i}
Allocate all tasks {T i ′ ,j ′ ,k ′ with i ′ = i or j ′ = j or k ′ = k} that are not yet processed to Pu and mark them processed has been done for Algorithm DynamicOuter in Section 3.3. The analysis of both processes is in fact rather similar, so that we will mostly state the corresponding lemmas, the proofs being similar to those presented in Section 3.

3.

In what follows, we will assume that N , the size of matrices A, B and C, is large and we will consider a continuous dynamic process whose behavior is expected to be close to the one of DynamicMatrix2Phases. In what follows, as in Section 3.3, we will concentrate on processor P k whose speed is s k and relative speed rs k = s k i s i . We will also denote by C = A × B the result of the matrix multiplication. Note that throughout this section, A i,k denotes the element of A on the ith row and jth column.

Let us assume that there exist 3 index sets I, J and K such that

• P k knows all elements A i,k , B k,j and Ci,j for any (i, j, k) ∈ I × J × K.

• I, J and K have size y.

In Algorithm DynamicMatrix2Phases, at each step, P k chooses to increase its knowledge by increasing y by l, which requires to receive (2y + 1)l elements of each A, B and C. As we did in Section 3.3, we will concentrate on x = y/N , and assuming that N is large, we will change the discrete process into a continuous process described by an ordinary differential equation depicting the evolution of expected values and we will rely on extensive simulations to assert that this approximation is valid.

In this context, let us consider that an elementary task T (i, j, k) consists in computing Ci,j ← Ci,j+A i,k B k,j . There are N 3 such tasks. In what follows, we will denote by g k (x) the fraction of elementary tasks that have not been computed yet at the instant when P k knows x 2 elements of A, B and C respectively, in the computational domain that does not include the tasks T (i, j, k) such that (i, j, k) ∈ I × J × K (this domain is equivalent to the "L"-shaped area for the outer product in Section 3.3). The following lemma enables to understand the dynamics of g k (all proofs are omitted because they are very similar to those of Section 3.3).

Lemma 7. g k (x) = (1 -x 3 ) α k , where α k = i =k s i s k .
Let us now denote by t k (x) the time step such that index sets I, J and K have size x. Then,

Lemma 8. t k (x) i si = 1 -N 2 (1 -(1 -x 3 ) α k +1 ).
Above equations well describe the dynamics of Dynamic-Matrix2Phases as long as it is possible to find elements of A, B and C that enable to compute enough unprocessed elementary tasks. On the other hand, as in the case of DynamicOuter2Phases, at the end, it is better to switch to another algorithm, where unprocessed elementary tasks T (i, j, k) are picked up randomly, what requires possibly to send all three values of A i,k , B k,j and Ci,j. In order to decide when to switch from one strategy to the other, let us introduce the additional parameter β.

As in the outer-product problem, a lower bound on the communication volume received by P k can be obtained by considering that each processor has a cube of tasks T i,j,k to compute, proportional to its relative speed. The edge-size of this cube is thus N 3 √ rs k . To compute all tasks in this cube, P k needs to receive a square of each matrix, that is 3N 2 rs 2/3 k . In order to determine when we should switch between Phase 1 and Phase 2, we can observe that if

x 3 k = βrs k - β 2 /2rs 2 k , then t k (x k ) i si = N 2 (1 -e -β (1 + o(rs k ))), so that at first order, t k (x k ) is independent of k. The instant t = N 2 i s i (1 -e -β
) is therefore chosen to switch between Phases 1 and 2.

As in the context of the outer product, we need to find the value of β that minimizes the volume of communications. If the switch occurs at time t = N 2 i s i (1 -e -β ), then • the volume of communications during Phase 1 is given by

3N 2 β 2/3 k rs 2/3 k -3N 2 β 5/3 k rs 5/3 k ,
• the volume of communications during Phase 2 is given by e

-β N 3 1 -β 2/3 k rs 5/3 k
, so that the total amount of communications with respect to the lower bound 3N 2 rs 2/3 k is given by

β 2/3 -β 5/3 k rs 5/3 k k rs 2/3 k + e -β N k rs 5/3 k 1 -β 2/3 k rs 5/3 k .

Simulation Results

We have conducted extensive simulations to compare the performance of the dynamic strategies with the previous analysis. Figure 9 presents the results for matrices of size 40x40 and Figure 10 presents the results for matrices of size 100x100. As in previous simulations, processor speeds are chosen uniformly at random in the interval [START_REF] Gautier | Kaapi: A thread scheduling runtime system for data flow computations on cluster of multi-processors[END_REF]100] and all amounts of communications have been normalized using the lower bound 3N 2 k rs 2/3 k on communications presented in the previous section.

As for the outer-product problem, we notice that dataaware strategies largely outperform simple strategies, and that DynamicMatrix2Phases is able to reduce the communication amount even more than DynamicMatrix. When the number of processors is large enough (i.e. in our simulation setting, p ≥ 50), our previous analysis is able to very accurately predict the performance of DynamicMatrix-2Phases.

Normalized communication amount 8 10

Value of β We also performed simulations of DynamicMatrix2Phases with varying values of β to check if the optimal value determined in the theoretical analysis actually minimizes the amount of communications. This is illustrated in Figure 11, for 100 processors, N/l = 40 and a fixed distribution of computing speeds. As for the outer product, we notice that the analysis accurately models the amount of communications of DynamicMatrix2Phases in the range of values of interest of β, and that the optimal value of β for the analysis (2.95) allows to obtain an amount of communications that is close to optimal. This corresponds to 94.7% of the tasks to be processed by the first phase of the algorithm. As for the outer product, we also notice that the value of β given by an analysis which is agnostic to processor speeds and assumes homogeneous speeds is very close to the optimal value (2.92 on this example).

CONCLUSION AND PERSPECTIVES

The contributions of this paper follow two directions. First, we have proposed randomized dynamic scheduling strategies for the outer product and the matrix multiplication kernels. We have proved that dynamic scheduling strategies that aim to place tasks on processors such that the induced amount of communications is as small as possible perform well. Second, we have been able to propose an Ordinary Differential Equation (ODE) whose solution describes very well the dynamics of the system. Even more important, we prove that the analysis of the dynamics of the ODE can be used in order to tune parameters and to inject some static knowledge which is useful to increase the efficiency of dynamic strategies. A lot remains to be done in this domain, that we consider as crucial given the practical and growing importance of dynamic runtime schedulers. First, it would be of interest to be able to provide analytical models for a larger class of dynamic schedulers even in the case of independent tasks, and to analyze their behavior also in dynamic environments (when the performance of the resources is unknown and varies over time). Then, it would be very useful to extend the analysis to applications involving both data and precedence dependencies. Extending this work to regular dense linear algebra kernels such as Cholesky or QR factorizations would be a promising first step in this direction.

ACKNOWLEDGEMENT

This work has been partially supported by the ANR SOL-HAR project, funded by the French Research Agency.

Figure 1 :

 1 Figure 1: Comparison of random and data-aware dynamic strategies, for vectors of size N/l = 100 blocks

Figure 2 :

 2 Figure 2: Communication amount of Dynamic-Outer2Phases and comparison to the other schedulers for different thresholds (for a given distribution of computing speeds with 20 processors and N/l = 100).

Figure 3 :

 3 Figure3: Illustration for the proof of Lemma 1. The top-left blue rectangle represents the data owned by the processor at time t k (x) (a permutation of the rows and columns has been applied to have it in the upper left corner). The new elements δx are depicted in red, as well as the corresponding available tasks. Note that some tasks (in black) corresponding to the combination of δx with the known elements have already been processed by other processors.

Figure 4 :

 4 Figure 4: Communication amounts of all outerproduct strategies for vectors of size N/l = 100 blocks ((N/l) 2 tasks).

Figure 5 :

 5 Figure 5: Communication amounts of all outerproduct strategies for vectors of size N/l = 1000 blocks ((N/l) 2 tasks).

Figure 6 :

 6 Figure 6: Communication amounts of Dynamic-Outer2Phases and its analysis for varying value of the β parameter which defines the threshold.

Figure 7 :

 7 Figure 7: Behavior of the heuristics for outer product for different values of heterogeneity (p = 20 processors and N/l = 100 blocks). For a given value h of heterogeneity, processor speeds are taken uniformly at random in the interval [100 -h, 100 + h].

Figure 8 :

 8 Figure 8: Behavior of the heuristics for outer product for different scenarios of heterogeneity (p = 20 processors and N/l = 100 blocks).

Figure 11 :

 11 Figure 11: Communication amount of Dynamic-Matrix2Phases and its analysis for varying value of the β parameter which defines the threshold.

Figure 9 :

 9 Figure 9: Communication amounts of all strategies for matrices of size N/l = 40 blocks (N 3 /l 3 = 64, 000 tasks).

Figure 10 :

 10 Figure 10: Communication amounts of all strategies for matrices of size N/l = 100 blocks (N 3 /l 3 = 1, 000, 000 tasks).


Youcef Mammeri
email: youcef.mammeri@u-picardie.fr

Yumeng Zhang

Comparison of solutions of Boussinesq systems

Keywords: Boussinesq systems, Comparison, Energy method. MS Codes. 35B60, 35Q51, 35Q53

We compare the solution of the generalized Boussinesq systems, for various values of a,b,c,d,

, that respective solutions of Boussinesq systems, starting from the same initial datum, remain close on a time interval inversely proportional to the wave amplitude.

Introduction

In this paper, we are interested in the two-way propagation of small amplitude long waves in shallow water. This phenomenon is described by the amplitude η and the velocity u of the waves, which satisfy the following generalized Boussinesq systems [START_REF] Bona | Boussinesq equations and other systems for smallamplitude long waves in nonlinear dispersive media I: Derivation and the linear theory[END_REF][START_REF] Bona | Boussinesq equations and other systems for smallamplitude long waves in nonlinear dispersive media II: Nonlinear theory[END_REF], for (x, t) ∈ R × R,

η t + u x + ε((ηu) x + au xxx -bη xxt ) = 0 u t + η x + ε(uu x + cη xxx -du xxt ) = 0.
Here 0 < ε < 1 denotes the quotient between the amplitude and the depth of water, and ε is assumed to be small. This coefficient is also proportional to the square of the quotient between the depth of water and the wavelength.

If we focus on the one-way propagation, the Boussinesq systems are reduced to the Korteweg-de Vries equation (KdV) [START_REF] Korteweg | On the change of form of long waves advancing in a rectangular canal and on a new type of long stationnary waves[END_REF] η t + η x + ε(α ηη x + β η xxx ) = 0, or the Benjamin-Bona-Mahony equation (BBM) [START_REF] Benjamin | Model equations for long waves in nonlinear dispersive systems[END_REF] η x + η t + ε(α ηη x -β η xxt ) = 0.

1 Recall that Bona, Pritchard and Scott [START_REF] Bona | A comparison of solutions of two model equations for long waves[END_REF] proved the following result.

Theorem 0.1 Let m ≥ 0 and f ∈ H m+5 (R). Let ε > 0 and η, ζ respectively the unique solution of KdV and BBM in C([0, ∞[; H m+5 (R)), starting from the same initial datum η| t=0 = ζ| t=0 = f . Then, there exists a constant M i > 0 with 0 ≤ i ≤ m such that for |t| ≤ ε -1 T , we have

||∂ i x η(., t) -∂ i x ζ(., t)|| L 2 (R) ≤ M i ε 2 |t|,
with T := Cm+5 f m+5 .

Concerning the Boussinesq systems, our result reads as follows.

Theorem 0.2 Let m ≥ 0 and (f, g) ∈ H m+5 (R) × H m+5 (R). Let (η 1 , u 1 ), (η 2 , u 2 ) be solution in C([-ε -1/2 T, ε -1/2 T ]; H m+5 (R)) 2 of Boussinesq system with (a 1 , b 1 , c 1 , d 1 ), respectively (a 2 , b 2 , c 2 , d 2 ), starting from the same initial datum (f, g) with

T := C m+5 (f, g) m+5 .
Then, for all 0 ≤ i ≤ m, there exists a constant M i > 0, depending only on (f, g) m+5 , such that for |t| ≤ ε -1/2 T , we have

∂ i x η 1 (•, t) u 1 (•, t) -∂ i x η 2 (•, t) u 2 (•, t) L 2 ≤ M i ε|t|.
This result reflects the fact that the solutions of Boussinesq systems with various coefficients are close with a maximum deviation of the order of ε 1/2 . This is not as good as the order of the neglected terms in the asymptotic expansion of the Euler equation. On the other hand, the obtained time of comparison is of order ε -1/2 , although the effects of dispersion and the effects of the nonlinearity have an order 1 influence. One can certainly hope that the solutions remain close longer, and improve the result, for example by choosing a suitable velocity, or a good combination of KdV solutions. In this way, Alazman et al. [START_REF] Alazman | Comparisons between the BBM equation and a Boussinesq system[END_REF] recently obtained a similar result and showed that the solution of the Boussinesq system with one-way propagation is well approximated by the solution of the KdV or BBM. Related papers concern the comparison of the solution of twodimensional Euler systems with the solution of Boussinesq systems or KdV equation [START_REF] Bona | Long wave approximations for water waves[END_REF][START_REF] Craig | An existence theory for water waves, and Boussinesq and Korteweg-de Vries scaling limits[END_REF][START_REF] Schneider | The long-wave limit for the water wave problem. I. The case of zero surface tension[END_REF][START_REF] Wright | Corrections to the KdV approximation for water waves[END_REF].

Our paper is organised as follows. In the first section, the local well-posedness of the Boussinesq systems is summed up. The second section deals with the comparison of the solution of the systems. Finally, we present some numerical illustrations.

Summary of existence theory

We consider the initial value problem, for (x, t) ∈ R × R,

η t + u x + ε ((ηu) x + au xxx -bη xxt ) = 0 (1.1)
u t + η x + ε (uu x + cη xxx -du xxt ) = 0 (1.2) η(x, 0) = f (x), u(x, 0) = g(x). (1.3)
The coefficients a, b, c, d satisfy one of these two conditions:

• b ≥ 0, d ≥ 0, a ≤ 0, c ≤ 0 • b ≥ 0, d ≥ 0, a = c > 0.
Recall that the Cauchy problem is locally well-posed [START_REF] Bona | Boussinesq equations and other systems for smallamplitude long waves in nonlinear dispersive media I: Derivation and the linear theory[END_REF][START_REF] Bona | Boussinesq equations and other systems for smallamplitude long waves in nonlinear dispersive media II: Nonlinear theory[END_REF]. In this work, we only need smooth data, and we refer to these papers for more precise assumptions on the regularity.

Theorem 1.1 Let m > 3/2 and (f, g) ∈ H m (R) × H m (R).
Then there exists C m > 0, depending only on m, such that for

T := C m (f, g) m ,
there exists a unique solution

(η, u) ∈ C([-ε -3/2 T, ε -3/2 T ]; H m (R))×C([-ε -3/2 T, ε -3/2 T ]; H m (R)) of the Cauchy problem (1.1)-(1.2)-(1.3). Moreover, let ( f , g) ∈ H m (R)×H m (R) be the initial datum of the unique solution ( η, u) and M > 0 such that (f, g) m ≤ M and ( f , g) m ≤ M.
Then there exists c m > 0 such that for

|t| ≤ ε -3/2 C m /M , η(•, t) u(•, t) - η(•, t) u(•, t) m ≤ c m f g - f g m .
Proof. The proof in [START_REF] Bona | Boussinesq equations and other systems for smallamplitude long waves in nonlinear dispersive media II: Nonlinear theory[END_REF] is actually done with ε = 1. The change of variables

x = ε -1/2 x, t = ε -1/2 t, η = ε -1 η, u = ε -1 u
provides the following initial value problem

η t + u x + ( η u) x + a u x x x -b η x x t = 0 u t + η x + u u x + c η x x x -d u x x t = 0 η( x, 0) = εf ( x), u( x, 0) = εg( x).
2 Comparison of solutions

The linear case

We consider the linearized problem around (0, 0)

η t + u x + ε(au xxx -bη xxt ) = 0 u t + η x + ε(cη xxx -du xxt ) = 0 η(x, 0) = f (x), u(x, 0) = g(x).
Let us compare the solution of the KdV-KdV system (a = c = 1, b = d = 0) and the solution of the BBM-BBM system (a = c = 0, b = d = 1). Assume that f, g ∈ S(R). The solution of KdV-KdV is given by the Fourier transform in space

η u (k, t) = e -ikA1(k)t f g (k), with A 1 (k) = 0 1 -εk 2 1 -εk 2 0 .
And the solution of BBM-BBM is

ζ v (k, t) = e -ikA2(k)t f g (k), with A 2 (k) = 0 1 1+εk 2 1 1+εk 2 0 . Since A 1 (k)A 2 (k) = A 2 (k)A 1 (k), it provides η u (k, t) - ζ v (k, t) = e -ikA1(k)t -e -ikA2(k)t f g (k) = e -ikA2(k)t e -ik(A1(k)-A2(k))t -1 f g (k). Thus η u (k, t) - ζ v (k, t) = e -ikA2(k)t e -ik(A1(k)-A2(k))t -1 f g (k) ≤ e -ik(A1(k)-A2(k))t -1 f g (k) .
However | exp(ix) -1| ≤ |x|, and we find

η u (k, t) - ζ v (k, t) ≤ -k(A 1 (k) -A 2 (k))t f g (k) .
Then, using 1 + εk 2 ≥ 1,

η u (., t) - ζ v (., t) ≤ ε 2 k 5 t 1 + εk 2 g f (k) ≤ ε 2 |k| 5 |t| g f (k) , and
η u (., t) - ζ v (., t) L ∞ ≤ 1 2π η u (., t) - ζ v (., t) L 1 ≤ ε 2 |t| R |k| 5 g f (k) dk.
The result can be written as follows.

Theorem 2.1 Let f and g be real valued functions with ∂ 5 x f and ∂ 5

x g integrables on R. Let (η, u) and (ζ, v) be solutions of the linear Cauchy problem associated with KdV-KdV and BBM-BBM respectively, starting from the same initial datum

η(•, 0) u(•, 0) = ζ(•, 0) v(•, 0) = f g .
Then there exists a constant C > 0, depending only on (f, g), such that, for all time t ∈ R,

η(•, t) u(•, t) - ζ(•, t) v(•, t) L ∞ ≤ Cε 2 |t|.

The main result

To simplify the reading, we only prove the result for the KdV-KdV and BBM-BBM system. The proof of the theorem 0.2 can be done in a similar way. We now consider the following nonlinear Cauchy problems. The first one associated with KdV-KdV:

η t + u x + ε ((ηu) x + u xxx ) = 0 (2.1)
u t + η x + ε(uu x + η xxx ) = 0 (2.2) η(x, 0) = f (x), u(x, 0) = g(x).
And the second one is called BBM-BBM:

ζ t + v x + ε((ζv) x -ζ xxt ) = 0 (2.3) v t + ζ x + ε(vv x -v xxt ) = 0 (2.4) ζ(x, 0) = f (x), v(x, 0) = g(x).
We prove that the respective solutions remain close as follows.

Theorem 2.2 Let m ≥ 5 and (f,

g) ∈ H m (R) × H m (R). Let (η, u) and (ζ, v) be solution in C([-ε -1/2 T, ε -1/2 T ]; H m (R))×C([-ε -1/2 T, ε -1/2 T ]; H m (R))
of KdV-KdV and BBM-BBM respectively, where

T := C m (f, g) m .
Then there exists a constant M > 0, depending only on

(f, g) m , such that for |t| ≤ ε -1/2 T , we have η(•, t) u(•, t) - ζ(•, t) v(•, t) L 2 ≤ M ε|t|.
Proof. To simplify the writings, we denote the L 2 -norm as ||.||. We assume that t > 0, negative time being dealt with similarly. Let us define

λ = η -ζ w = u -v.
Then (λ, w) satisfies:

λ t + w x -ελ xxt = -ε ((λw) x + (λu + ζw) x + u xxx + η xxt ) (2.5)
w t + λ x -εw xxt = -ε (ww x + (vw) x + η xxx + u xxt ) (2.6) λ(x, 0) = 0, w(x, 0) = 0.
Equations (2.5) and (2.6) are multiplied by λ and w respectively, and are integrated over space.

+∞ -∞ λ(λ t + w x -ελ xxt ) + w(w t + λ x -εw xxt ) dx = -ε +∞ -∞ λ(u xxx + η xxt ) + w(η xxx + u xxt ) -ε +∞ -∞ λ((λw) x + (λu + ζw) x ) + w(ww x + (vw) x ) dx.
(2.7) We first simplify the left-hand side of (2.7). Since λ, w ∈ H m (R), we know that, for 0 ≤ i ≤ m, ∂ i

x λ and ∂ i x w go to zero at infinity. Thus

+∞ -∞ (λw x + wλ x )dx = +∞ -∞ (λw) x dx = [λw] +∞ -∞ = 0,
and an integration by parts gives

+∞ -∞ λ(λ t -ελ xxt ) + w(w t -εw xxt )dx = 1 2 d dt +∞ -∞ λ 2 + w 2 + ελ 2 x + εw 2 x dx.
Integration by parts are applied to the right-hand side of (2.7), and it gets:

+∞ -∞ λ((λw) x + (λu + ζw) x )dx = - +∞ -∞ λ x (λw + λu + ζw) dx = - +∞ -∞ λ 2 2 x (w + u) dx - +∞ -∞ λ x ζw dx = +∞ -∞ λ 2 2 (w x + u x ) dx - +∞ -∞ λ x ζw dx. Similarly +∞ -∞ w(vw) x dx = - +∞ -∞ w x vw dx = - +∞ -∞ w 2 2 x v dx = +∞ -∞ w 2 2 v x dx,
and

+∞ -∞ w 2 w x dx = +∞ -∞ w 3 3 x dx = w 3 3 +∞ -∞ = 0.
Finally, the equality (2.7) becomes

1 2 d dt +∞ -∞ (λ 2 + w 2 + ελ 2 x + εw 2 x ) dx = -ε +∞ -∞ λ 2 2 (w x + u x ) + w 2 2 v x + λ x ζw dx -ε +∞ -∞ λu xxx + λη xxt + wη xxx + wu xxt dx.
However λ(x, 0) = 0 and w(x, 0) = 0, we deduce

+∞ -∞ (λ 2 + w 2 + ελ 2 x + εw 2 x ) dx ≤ 2ε t 0 +∞ -∞ λ 2 2 (w x + u x ) + w 2 2 v x + |λ x ζw| + 2ε t 0 +∞ -∞ |λu xxx | + |λη xxτ | + |wη xxx | + |wu xxτ |dxdτ.
(2.8)

According to theorem 1.1, we have for t

∈ [0, ε -1/2 T ], h = η, u, ζ, or v, and 0 ≤ i ≤ m, ∂ i x h ≤ h m ≤ C 0 = C 0 (f, g), moreover, for 0 ≤ i ≤ m -1 and t ∈ [0, ε -1/2 T ], ∂ i x h ∞ ≤ √ 2( ∂ i x h ∂ i+1 x h ) 1 2 ≤ √ 2C 0 .
Let us find an upper bound of the inequality (2.8). On one hand, we have

t 0 +∞ -∞ λ 2 2 (w x + u x ) dxdτ ≤ 1 2 sup τ ∈[0,ε -1/2 T ] ( u x (τ ) ∞ + w x (τ ) ∞ ) ≤ 1 2 sup τ ∈[0,ε -1/2 T ] (2 u x (τ ) ∞ + v x (τ ) ∞ ) t 0 λ 2 dτ ≤ 3 √ 2 2 C 0 t 0 λ 2 dτ, t 0 +∞ -∞ w 2 2 v x dxdτ ≤ 1 2 sup τ ∈[0,ε -1/2 T ] v x (τ ) ∞ t 0 w 2 dτ ≤ √ 2 2 C 0 t 0 w 2 dτ,
and, using the Cauchy-Schwarz inequality,

t 0 +∞ -∞ |λ x ζw|dxdτ ≤ sup τ ∈[0,ε -1/2 T ] ζ x ∞ t 0 λ x w dτ ≤ √ 2C 0 t 0 λ x w dτ.
On the other hand, the Cauchy-Schwarz inequality provides

t 0 +∞ -∞ |λu xxx |dxdτ ≤ t 0 λ u xxx dτ ≤ C 0 t 0 λ dτ t 0 +∞ -∞ |wη xxx |dxdτ ≤ t 0 w η xxx dτ ≤ C 0 t 0 w dτ.
However,

η xxt = -∂ 2 x (u x + ε((ηu) x + u xxx )) = -(u xxx + 3εη xx u x + 3εη x u xx + εηu xxx + εη xxx u + εu xxxxx ), thus, for t ∈ [0, ε -1/2 T ], η xxt ≤ u xxx + 3ε η xx u x ∞ + 3ε η x ∞ u xx + ε η ∞ u xxx + ε η xxx u ∞ + ε u xxxxx ≤ (1 + ε)C 0 + ε8 √ 2C 2 0 .
Thanks to the Cauchy-Schwarz inequality, we obtain

t 0 +∞ -∞ |λη xxτ |dxdτ ≤ t 0 λ η xxτ dτ ≤ (1 + ε)C 0 + ε8 √ 2C 2 0 t 0 λ dτ.
In the same way, we find for t ∈ [0, ε -1/2 T ],

u xxt = -∂ 2 x (η x + ε(uu x + η xxx )) ≤ ( η xxx + 3ε u x ∞ u xx + ε u ∞ u xxx + ε η xxxxx ) ≤ (1 + ε)C 0 + ε4 √ 2C 2 0 , and then t 0 +∞ -∞ |wu xxτ |dxdτ ≤ t 0 w u xxτ dτ ≤ (1 + ε)C 0 + ε4 √ 2C 2 0 t 0 w dτ.
The inequality (2.8) is finally rewritten

+∞ -∞ (λ 2 + w 2 + ελ 2 x + εw 2 x ) dx ≤ 2ε 3 √ 2 2 C 0 t 0 λ 2 dτ + √ 2 2 C 0 t 0 w 2 dτ + √ 2C 0 t 0 λ x w dτ + 2ε (2 + ε)C 0 + ε8 √ 2C 2 0 t 0 λ dτ + (2 + ε)C 0 + ε4 √ 2C 2 0 t 0 w dτ .
We define

A(t) := +∞ -∞ (λ 2 + w 2 + ελ 2 x + εw 2 x )dx 1/2 , thus λ 2 + w 2 ≤ A 2 and λ x 2 + w x 2 ≤ ε -1 A 2 ,
and the inequality (2.8) becomes, using 2ab ≤ a 2 + b 2 ,

A 2 (t) ≤ 2C 1 ε t 0 A(τ )dτ + C 1 ε 1/2 t 0 A 2 (τ )dτ ≤ 3C 1 ε t 0 A(τ )dτ + C 1 t 0 A 3 (τ )dτ
The following Gronwall's lemma is used [START_REF] Albert | Comparisons between model equations for long waves[END_REF].

Lemma 2.3 Let α, β > 0 and ρ > 1. Define T := β -1/ρ α (1-ρ)/ρ +∞ 0 (1 + z ρ ) -1 dz
Then there exists a constant C = C(ρ) > 0, independent of α and β such that if A is a non-negative, continuous function defined on [0, T ] with A(0) = 0 and for all t ∈ [0, T ]

A 2 (t) ≤ t 0 αA(τ ) + βA ρ+1 (τ )dτ,
then for all t ∈ [0, T ], A(t) ≤ Cαt.

Even if we choose T smaller, we conclude that, for 0

≤ t ≤ ε -1/2 T , η(•, t) u(•, t) - ζ(•, t) v(•, t) L 2 ≤ Cεt.
We can generalize this result to the space derivatives as follows.

Theorem 2.4 Let m ≥ 0 and (f,

g) ∈ H m+5 (R) × H m+5 (R). Let (η, u) and (ζ, v) be solution in C([-ε -1/2 T, ε -1/2 T ]; H m+5 (R)) × C([-ε -1/2 T, ε -1/2 T ]; H m+5 (R))
of KdV-KdV and BBM-BBM respectively, where

T := C m+5 (f, g) m+5 .
Then, for all 0 ≤ i ≤ m, there exists a constant M i > 0, depending only on (f, g) m+5 , such that for |t| ≤ ε -1/2 T , we have

∂ i x η(•, t) u(•, t) -∂ i x ζ(•, t) v(•, t) L 2 ≤ M i ε|t|.
Proof. The proof is done by induction on the derivative order. The preceding theorem gives the first step i = 0, and we assume that, for 0 ≤ i ≤ n -1 and 0 ≤ t ≤ ε -1/2 T , there exists

M i = M i ( (f, g) (m+5)×(m+5) ) > 0 such that ∂ i x λ L 2 + ∂ i x w L 2 = ∂ i x η(•, t) u(•, t) -∂ i x ζ(•, t) v(•, t) L 2 ≤ M i εt.
The result is now shown for the step n. Multiplying equations (2.5) and (2.6) by ∂ 2n

x λ and ∂ 2n x w respectively, and integrating over space:

+∞ -∞ ∂ 2n x λ(λ t + w x -ελ xxt ) + ∂ 2n x w(w t + λ x -εw xxt ) dx = -ε +∞ -∞ ∂ 2n x λ(u xxx + η xxt ) + ∂ 2n x w(η xxx + u xxt ) -ε +∞ -∞ ∂ 2n x λ((λw) x + (λu + ζw) x ) + ∂ 2n x w(ww x + (vw) x ) dx.
(2.9) Concerning the left-hand side of (2.9), integration by parts provides

+∞ -∞ (∂ 2n x λ w x + ∂ 2n x w λ x )dx = (-1) n +∞ -∞ (∂ n x λ∂ n x w) x dx = 0,
and

+∞ -∞ ∂ 2n x λ(λ t -ελ xxt )+∂ 2n x w(w t -εw xxt )dx = (-1) n 2 d dt +∞ -∞ (∂ n x λ) 2 +(∂ n x w) 2 +ε(∂ n+1 x λ) 2 +ε(∂ n+1 x w 2 )dx.
The right-hand side of (2.9) is rewritten as follows.

-

(-1) n ε +∞ -∞ ∂ n x λ ∂ n+1 x (λw + λu + ζw) + ∂ n x w ∂ n x (ww x ) + ∂ n x w ∂ n+1 x (vw) dx -(-1) n ε +∞ -∞ ∂ n x λ ∂ n+3 x u + ∂ n x λ ∂ n+2 x η t + ∂ n x w ∂ n+3 x η + ∂ n x w ∂ n+2 x u t dx.
Let us define

B 2 (t) := +∞ -∞ (∂ n x λ) 2 + (∂ n x w) 2 + ε(∂ n+1 x λ) 2 + ε(∂ n+1 x w) 2 dx.
Since λ(•, 0) = 0 and w(•, 0) = 0, the equation (2.9) becomes

B 2 (t) = -2ε t 0 +∞ -∞ ∂ n x λ ∂ n+1 x (λw + λu + ζw) + ∂ n x w ∂ n x (ww x ) + ∂ n x w ∂ n+1 x (vw) dxdτ -2ε t 0 +∞ -∞ ∂ n x λ ∂ n+3 x u + ∂ n x λ ∂ n+2 x η τ + ∂ n x w ∂ n+3 x η + ∂ n x w ∂ n+2 x u τ dxdτ.
Thanks to the Leibniz formula, we find

t 0 +∞ -∞ ∂ n x λ ∂ n+1 x (λu)dxdτ = t 0 +∞ -∞ ∂ n x λ n+1 i=0 n + 1 i ∂ n+1-i x λ ∂ i x u dxdτ = t 0 +∞ -∞ ∂ n x λ ∂ n+1 x λ u + (n + 1)(∂ n x λ) 2 u x dxdτ + t 0 +∞ -∞ n+1 i=2 n + 1 i ∂ n x λ ∂ n+1-i x λ ∂ i x udxdτ. However t 0 R ∂ n x λ ∂ n+1 x λ udxdτ = t 0 R 1 2 ∂ x (∂ n x λ) 2 udxdτ = - t 0 R 1 2 (∂ n x λ) 2 u x dxdτ, then we have t 0 R ∂ n x λ ∂ n+1 x (λu)dxdτ = t 0 R (n + 1 2 )(∂ n x λ) 2 u x + n+1 i=2 n + 1 i ∂ n x λ ∂ n+1-i x λ ∂ i x u dxdτ.
In the same manner, it gets

t 0 R ∂ n x w ∂ n+1 x (vw)dxdτ = t 0 R (n + 1 2 )(∂ n x w) 2 v x + n+1 i=2 n + 1 i ∂ n x w ∂ n+1-i x w ∂ i x v dxdτ,
and t 0 R ∂ n x w ∂ n x (ww x )dxdτ = t 0 R (2n - 1 2 )(∂ n x w) 2 w x + n-1 i=2 n i ∂ n x w ∂ i x w ∂ n+1-i x w dxdτ.
Thus,

B 2 (t) = -2ε t 0 R (n + 1 2 )(∂ n x λ) 2 (w x + u x ) + (n + 1 2 )(∂ n x w) 2 v x + (2n - 1 2 )(∂ n x w) 2 w x + ∂ n x λ ∂ n+1 x ζ∂ n x w + ∂ n x λ ∂ n+1 x w∂ n x ζ + n-1 i=2 n i ∂ n x w ∂ i x w ∂ n+1-i x w + n+1 i=2 n + 1 i (∂ n x λ ∂ n+1-i x λ (∂ i x w + ∂ i x u) + ∂ n x w∂ n+1-i x w ∂ i x v + ∂ n x λ∂ n+1-i x w ∂ i x ζ) +∂ n x λ ∂ n+3 x u + ∂ n x λ ∂ n+2 x η τ + ∂ n x w ∂ n+3 x η + ∂ n x w ∂ n+2 x u τ dxdτ. (2.10) Recall that, for t ∈ [0, ε -1/2 T ], h = η, u, ζ, v and 0 ≤ i ≤ m + 4, we have ∂ i x h ∞ ≤ √ 2( ∂ i x h ∂ i+1 x h ) 1 2 ≤ √ 2C 0 .
We deduce

t 0 R (n + 1 2 )(∂ n x λ) 2 (w x + u x ) dxdτ ≤ (n + 1 2 ) t 0 ∂ n x λ 2 ( w x ∞ + u x ∞ )dτ ≤ (n + 1 2 )2 √ 2C 0 t 0 ∂ n x λ 2 dτ,
(2.20 -a)

similarly t 0 R (n + 1 2 )(∂ n x w) 2 v x dxdτ ≤ (n + 1 2 ) √ 2C 0 t 0 ∂ n x w 2 dτ, (2.20 -b) t 0 R (2n - 1 2 )(∂ n x w) 2 w x dxdτ ≤ (2n - 1 2 )2 √ 2C 0 t 0 ∂ n x w 2 dτ, (2.20 -c) and t 0 R |∂ n x λ ∂ n+1 x ζ∂ n x w+∂ n x λ ∂ n+1 x w∂ n x ζ|dxdτ ≤ √ 2C 0 t 0 ∂ n x λ ∂ n x w dτ + t 0 ∂ n x λ ∂ n+1 x w dτ .
Using the Cauchy-Schwarz inequality, we obtain

t 0 R n+1 i=2 n + 1 i ∂ n x λ∂ n+1-i x λ∂ i x u dxdτ ≤ n+1 i=2 n + 1 i t 0 ∂ n x λ ∂ n+1-i x λ ∂ i x u ∞ dτ,
and the inductive assumption provides, for 0

≤ t ≤ ε -1/2 T , t 0 R n+1 i=2 n + 1 i ∂ n x λ∂ n+1-i x λ∂ i x u dxdτ ≤ n+1 i=2 n + 1 i M n+1-i t 0 ετ ∂ n x λ dτ ≤ n+1 i=2 n + 1 i M n+1-i ε 1/2 T 1/2 t 0 ∂ n x λ dτ.

Same computations give

t 0 R n+1 i=2 n + 1 i |(∂ n x λ ∂ n+1-i x λ (∂ i x w + ∂ i x u) + ∂ n x w∂ n+1-i x w ∂ i x v + ∂ n x λ∂ n+1-i x w ∂ i x ζ)| ≤ 3 n+1 i=2 n + 1 i M n+1-i ε 1/2 T 1/2 t 0 ∂ n x λ + ∂ n x w dτ. and t 0 R n-1 i=2 n i ∂ n x w∂ i x w∂ n+1-i x w dxdτ ≤ n-1 i=2 n i M i ε 1/2 T 1/2 t 0 ∂ n x w dτ, (2.20 -f )
On the other hand, from the Cauchy-Schwarz inequality, it gets

t 0 R ∂ n x λ∂ n+3 x u dxdτ ≤ t 0 ∂ n x λ ∂ n+3 x u dτ ≤ √ 2C 0 t 0 ∂ n x λ dτ, t 0 R ∂ n x w∂ n+3 x η dxdτ ≤ t 0 ∂ n x w ∂ n+3 x η dτ ≤ √ 2C 0 t 0 ∂ n x w dτ. Since ∂ n+2 x η t = -∂ n+2 x (u x + ε(ηu) x + εu xxx ) = -(∂ n+3 x u + ε n+3 i=0 n + 3 i ∂ i x η ∂ n+3-i x u + ε∂ n+5 x u), it implies t 0 R ∂ n x λ∂ n+2 x η t dxdτ ≤ t 0 ∂ n x λ ∂ n+2 x η t dτ ≤ (1 + ε)C 0 + √ 2C 2 0 ε n+3 i=0 n + 3 i t 0 ∂ n x λ dτ.
In the same way, we have

t 0 R ε|∂ n x w∂ n+2 x u t |dxdτ ≤ t 0 ε ∂ n x w ∂ n+2 x u t dτ ≤ (1 + ε)C 0 + √ 2C 2 0 ε n+2 i=0 n + 2 i t 0 ∂ n x w dτ, (2.20 -k) because ∂ n+2 x u t = -∂ n+2 x (η x + εuu x + εη xxx ) = -(∂ n+3 x η + ε n+2 i=0 n + 2 i ∂ i x u ∂ n+3-i x u + ε∂ n+5 x η).
For t ∈ [0, ε -1/2 T ], the function B(t) being such that

∂ n x λ 2 + ∂ n x w 2 ≤ B 2 (t), and ∂ n+1 x λ 2 + ∂ n+1 x w 2 ≤ ε -1 B 2 (t),
the equation (2.10) finally becomes

B 2 (t) ≤ 2C n ε t 0 B(τ )dτ + C n ε 1/2 t 0 B 2 (τ )dτ ≤ 3C n ε t 0 B(τ )dτ + C n t 0 B 3 (τ )dτ.
We conclude using the Gronwall's lemma as for the initial step.

Corollary 2.5 Let m ≥ 1 and (f, g) ∈ H m+5 (R) × H m+5 (R). Let (η, u) and (ζ, v) be solution in C([-ε -1/2 T, ε -1/2 T ]; H m+5 (R)) × C([-ε -1/2 T, ε -1/2 T ]; H m+5 (R))
of KdV-KdV and BBM-BBM respectively, where

T := C m+5 (f, g) m+5 .
Then, for all 0 ≤ i ≤ m -1, there exists a constant N i > 0, depending only on (f, g) m+5 , such that for |t| ≤ ε -1/2 T , we have

∂ i x η(•, t) u(•, t) -∂ i x ζ(•, t) v(•, t) L ∞ ≤ N i ε|t|. Proof. Using ∂ i x (η -ζ) ∞ ≤ √ 2( ∂ i x (η -ζ) ∂ i+1 x (η -ζ) ) 1/2
, the preceding theorem is applied with 0 ≤ i ≤ m and |t| ≤ ε -1/2 T , to give

∂ i x (η -ζ) ∞ ≤ √ 2(M i ε|t|M i+1 ε|t|) 1/2 = 2M i M i+1 ε|t|, and
∂ i x (u -v) ∞ ≤ 2M i M i+1 ε|t|. It is enough to choose N i = 2 2M i M i+1 .
Remark 2.6 Concerning the generalized Boussinesq systems, the difference between the respective solutions can be written as follows:

λ t + w x + ε (a 2 w xxx -b 2 λ xxt ) = -ε ((λw) x + (λu + ζw) x + (a 1 -a 2 )u xxx + (b 2 -b 1 )η xxt ) w t + λ x + ε (c 2 λ xxx -d 2 w xxt ) = -ε (ww x + (vw) x + (c 2 -c 1 )η xxx + (d 2 -d 1 )u xxt ) λ(x, 0) = 0, w(x, 0) = 0, where λ = η 1 -η 2 , w = u 1 -u 2 .
The rest of the proof is similar to the previous cases.

Numerical comparison

The aim of this section is not to perform a complete numerical study, but rather to present some relevant simulations and to explore if the result is optimal and if the time comparison is longer than the theoretical one. Therefore we have an overview of the way in which solutions of each system evolve and differ. To simulate the system, a Crank-Nicolson scheme is used to discretize the time-derivative, and the fast Fourier transform is applied for the space derivatives. The obtained non-linear problem is solved using fixed point iterations.

We consider a bounded domain [-L, L], L > 0 a large fixed value. We denote N x > 0 the number of Fourier modes, ∆t > 0 the time step, and for n ∈ N, ηn , resp. ûn , is the approximation of η, resp. û, at time n∆t . Algorithm: By denoting by m ≥ 1 a maximal number of the Picard iterations, the algorithm is described as follows:

• Set η 0,m := ... := η 0,1 := η 0,0 := η 0 and u 0,m := ... := u 0,1 := u 0,0 := u 0 .

• For n = 0, 1, ..., compute:

• For r = 0, 1, ..., m -1,

I + ∆t 2 A(ξ) η n+1,r+1 u n+1,r+1 = I - ∆t 2 A(ξ) η n u n (3.1) - iεξ∆t 4
1 1+εbξ 2 (η n+1,r + η n )(u n+1,r + u n ) 1 2(1+εdξ 2 ) (u n+1,r + u n ) 2 , where, ξ = k π L , -Nx 2 ≤ k ≤ Nx 2 -1, I = 1 0 0 1 ; A(ξ) = iξ 0 1-εaξ 2 ) 1+εdξ 2 1-εcξ 2 1+εdξ 2 0 .
The iterations are stopped in one of the two following cases:

• when ηn+1,r+1-ηn+1,r l 2 ηn+1,0 l 2 + un+1,r+1-un+1,r l 2 un+1,0 l 2 ≤ τ , with τ > 0 a fixed tolerance.

• or when r = m -1. We mention that the step ∆t can be reduced in this case, in order to improve the previous relative error.

We then set η n+1 := η n+1,r+1 and u n+1 := u n+1,r+1 .

Remark 3.1 We choose a suitable step ∆t > 0 as follows:

∆t < ∆x πε(2 + ε(a + c)(π/∆x) 2 )( η n l ∞ + u n l ∞ ) . (3.2)
Indeed, from (3.1),

I + ∆t 2 A(ξ) η n+1 -η n+1,r+1 u n+1 -u n+1,r+1 = - iεξ∆t 4   1 1+εbξ 2 (η n+1 + η n )(u n+1 + u n ) -(η n+1,r + η n )(u n+1,r ) 1 2(1+εdξ 2 ) (u n+1 + u n ) 2 -(u n+1,r + u n ) 2   .
The nonlinear term can be written as

(ηn+1-ηn+1,r)(un+1+un+1,r) 2 + (ηn+1+ηn+1,r)(un+1-un+1,r) 2 + u n (η n+1 -η n+1,r ) + η n (u n+1 -u n+1,r ) (u n+1 -u n+1,r )(u n+1 + 2u n + u n+1,r )
and a linearization, near (η n , u n ), in the right-hand side of this relation leads to the approximation:

I + ∆t 2 A(ξ) η n+1 -η n+1,r+1 u n+1 -u n+1,r+1 ≈ - iεξ∆t 4   1 1+εbξ 2 u n (η n+1 -η n+1,r ) + η n (u n+1 -u n+1,r ) 1 2(1+εdξ 2 ) u n (u n+1 -u n+1,r )   .
The matrix

I + ∆t 2 A(ξ) is invertible if, for all ξ, ∆t < 2(1 + εbξ 2 )(1 + εdξ 2 ) |ξ(1 -εaξ 2 )(1 -εcξ 2 )| ,
and it follows according to the Parseval's formula,

η n+1,r+1 -η n+1 l 2 + u n+1,r+1 -u n+1 l 2 ε|ξ|∆t 4 1 1 + εbξ 2 1 + |1 -εcξ 2 | 1 + εdξ 2 ( u n l ∞ η n+1,r+1 -η n+1 l 2 + η n l ∞ u n+1,r+1 -u n+1 l 2 ) + 1 1 + εdξ 2 1 + |1 -εaξ 2 | 1 + εbξ 2 ( u n l ∞ u n+1,r+1 -u n+1 l 2 ) . However ξ = k π L , with -Nx 2 ≤ k ≤ Nx 2 -1, and N x = 2L ∆x , it follows η n+1,r+1 -η n+1 l 2 + u n+1,r+1 -u n+1 l 2 επ∆t 4∆x 2 + ε(a + c)(π/∆x) 2 ( η n l ∞ + u n l ∞ ) ( η n+1,r -η n+1 l 2 + u n+1,r+1 -u n+1 l 2 ).
Remark 3.2 When b, d > 0 and since |ξ|/(1 + ξ 2 ) ≤ 1, a better choice for the time step ∆t > 0 can be done as follows:

∆t < 1 ε( η n l ∞ + u n l ∞ )
.

Simulations are performed with such a time step ∆t, and with L = 100, N x = 2 10 . Similar results were obtained with smaller space steps. We start from the initial datum:

η 0 (x) = 15 4 -2 + cosh(3 2/5x) sech 4 3 √ 10 x , u 0 (x) = 3wsech 2 3 √ 10
x , with w = ±5/2. This initial datum provides a solitary wave solution of the BBM-BBM system with ε = 1, a = c = 0, b = d = 1/6, given by η(x, t) = η 0 (x -wt) and u(x, t) = u 0 (x -wt). To ensure the convergence of the scheme, the following conservation law is tested :

∀t ∈ R, +∞ -∞ η 2 (x, t)+u 2 (x, t)+εb∂ x η 2 (x, t)+εd∂ x u 2 (x, t)dx = +∞ -∞ η 2 0 (x)+u 2 0 (x)+εb∂ x η 2 0 (x)+εd∂ x u 2 0 (x)dx.
Figure 1 shows the evolution over time of the solution, the L ∞ -norm, the conservation law above, and the error between the approximate solution and the solitary wave. On the top left, the approximated amplitude η at time t = 0, 10, 50. On the top right, the approximated velocity u at time t = 0, 10, 50. On the bottom left, the H 1 and L ∞ -norms. On the bottom right, the error between the approximated solution and the solitary wave.

Based on these results, the numerical scheme appears to be relevant for the simulations. The approximate solution remains close to the exact soliton, and the conservation laws are well preserved.

We represent in Figure 2 the evolution of the solution starting from a localized wave defined as follows η 0 (x) = u 0 (x) = α exp(-x 2 /2 2 ). The solution of the KdV-KdV system disperses faster, the solution of the BBM-BBM system being the slowest. The difference between the solutions increases rapidly, the slowdown is due to the decay of solutions. Between t = 1/ √ ε, and t = 1/ε, the difference has doubled. We also noticed that conservation laws are well preserved.

We inspect in Figure 3 the influence of the parameter ε and of the amplitude of initial data. We find that the difference between solution increases with ε. Nevertheless, the growth is quickly done during the first iterations. On the other hand, the constants depending on the norm of the initial data, it appears that the difference between the solutions is even greater than the amplitude of the solution is great and the dispersion is slow.

To sum up, when α and ε are sufficiently small, the solutions behave like the linear case. It is due to the dispersion, and thus the smallness of the solutions. Then the solutions remain close and the result can be improved. If α and ε become large, the solutions differ quickly. Nevertheless, it is possible to improve the result by selecting special solutions, e.g. by choosing one-way propagation solution [START_REF] Alazman | Comparisons between the BBM equation and a Boussinesq system[END_REF][START_REF] Bona | A comparison of solutions of two model equations for long waves[END_REF].

Figure 1 :

 1 Figure 1: Result for ε = 1, a = c = 0, b = d = 1.On the top left, the approximated amplitude η at time t = 0, 10, 50. On the top right, the approximated velocity u at time t = 0, 10, 50. On the bottom left, the H 1 and L ∞ -norms. On the bottom right, the error between the approximated solution and the solitary wave.

Figure 2 :

 2 Figure 2: Result for ε = 0.1, α = 1. On the top, the approximated amplitude η at time t = 1/ √ ε, 1/ε, 1/ε 2 . At the center, the approximated velocity u at time t = 1/ √ ε, 1/ε, 1/ε 2 . On the bottom left, the conservation laws. On the bottom right, the difference between the approximated solution of the Boussinesq system (a = b = c = d = 1), the KdV-KdV system (a = c = 1, b = d = 0) and the BBM-BBM system (a = c = 0, b = d = 1).

Figure 3 :

 3 Figure 3: Difference between the approximated solution of the Boussinesq system (a = b = c = d = 1), the KdV-KdV system (a = c = 1, b = d = 0) and the BBM-BBM system (a = c = 0, b = d = 1). At left, result for ε = 0.5, 0.1, 0.01 and α = 1. At right, result for ε = 0.1 and α = 1, 2, 3.