Anne Etien

Vincent Aranega

Xavier Blanc
email: xavier.blanc@labri.fr

Richard F Paige
email: richard.paige@york.ac.uk

Chaining Model Transformations

Model transformation is one of the key practices of Model-Driven Engineering. Building very large model transformations may benefit from the construction of small transformations, in order to manage complexity and enhance reusability, maintainability and modularity. The decomposition of transformations into smaller ones raises the issue of assuring the validity of a composition: if two or more transformations are chained together, are the results of executing the chain the expected ones?

This paper addresses the challenge of determining if two transformations are conflicting. Transformations can conflict in numerous ways, e.g., in terms of preconditions, postconditions, behaviour of individual rules. In this paper, we demonstrate a strong notion of conflict, via commutativity: two transformations do not conflict if they can be chained in either order, and in doing so produce identical results. We propose an approach to detecting such potential conflicts based on static analysis, exploiting an intermediate representation of transformations independent of any concrete language.

INTRODUCTION

Model Driven Engineering (MDE) advocates the principle of separation of concerns, through the extensive use of models in all steps of the software development cycle [START_REF] Selic | The pragmatics of model-driven development[END_REF]. A number of approaches [START_REF] Pilgrim | Constructing and visualizing transformation chains[END_REF][START_REF] Vanhooff | A framework for transformation chain development processes[END_REF] recommend to adopt this principle for designing and implementing transformations. Thus, it is desirable to decompose transformations, much as engineers decompose other artefacts (like architectures and object-oriented designs). If large transformations are suitably decomposed, their scalability, maintainability and reusability may be improved [START_REF] Pilgrim | Constructing and visualizing transformation chains[END_REF][START_REF] Wagelaar | Module superimposition: a composition technique for rule-based model transformation languages[END_REF], and it thereafter becomes tractable for transformations to be systematically engineered [START_REF] Cordy | Eating our own dog food: Dsls for generative and transformational engineering[END_REF].

The decomposition of transformations into smaller ones raises a new issue: assuring the validity of the composition -i.e. the chaining -of these transformations. Understanding what constitutes a valid chain of transformations, and providing support for checking validity, is the focus of this paper.

This paper deals with transformations where the intersection of the input and output metamodels is not empty. In [START_REF] Etien | Combining independent model transformations[END_REF], we proposed an analysis based on the metamodels involved in the transformations to identify when transformations can be chained. These constraints effectively allow us to prune some chains that are invalid. However, this analysis had a limitation: in some situations where two transformations f and g were analysed, it was determined that these transformations could be chained in either order (e.g. f • g or g • f ), though each ordering produced different models. A finer grained analysis is thus required. This paper addresses the challenge of identifying conflicts between transformations and thus determining if any two arbitrary transformations are commutative, i.e. if they can be composed in either order and produce the same results no matter what the input model. The main contribution of the paper is to propose a static analysis of the transformations that relies on the following intuition: a transformation reading some instances of a metaclass is in conflict with another transformation creating, deleting or modifying instances of the same metaclass.

The paper is organized as follows. Section 2 presents a motivating example. Section 3 provides a set of formal definitions, in Alloy, that allow us to define what constitutes a valid chain of transformations, and thereafter derive associated chaining constraints. Section 4 details how the checking of the chaining constraints can be implemented. Section 5 illustrates the advantages of the proposed static analysis. Section 6 describes related work on transformation composition and reusability. Section 7 draws some conclusions and additional perspectives.

MOTIVATING EXAMPLE

Consider a scenario where we want to automatically generate an information system from a UML model. The generated system should satisfy established practices like adopting the 3-tier Model View Controller architecture. We have some reusable transformations available, including: a transformation T1 applying the MVC framework and a transformation T2 introducing an identifier to each class. Concretely, transformation T1 handles UML concepts and introduces View and Control and their relationships with Classes from UML in the output model, while transformation T2 adds an identifier property to each Class. Building a valid chain from these two transformations is not so simple and the task becomes harder when the number of transformations increases.

Intuitively, two transformations can be composed in sequence if the concepts used by the second transformation are not deleted by the first one. Applying this constraint in either order establishes if the two transformations can be inverted. These chaining constraints have been identified and formally defined in [START_REF] Etien | Combining independent model transformations[END_REF]. Thus, checking the combination in both orders, according to this approach, is equivalent to showing that the difference between, the intersection of the input metamodels of both transformations and the intersection of the output metamodels, is empty. More formally, if (InT 1 ∩ InT 2 ) \ (OutT 1 ∩ OutT 2 ) = ∅ then T1 and T2 can be chained in both order (with InT 1 (resp. InT 2 ) the input metamodel of T1 (resp. T2) and OutT 1 (resp. OutT 2 ) its output one). On the motivating example, we lead to the conclusion that T1 and T2 can be chained in both orders; however, the chaining constraints do not take into account whether the transformations are conflicting and thus commutative i.e. that the chaining T1 and T2 or T2 and T1 leads to the same results. In our scenario, it is clear that T1 and T2 are not commutative. In applying T1 then T2, all the classes (including the Control ) will have an ID, whereas in the opposite case, only the classes initially designed. The two transformations are conflicting. The resulting system will be different. The two orders are valid, but only one leads to the expected result. The purpose of this paper is to provide constraints to automatically prevent the chain designer that a verification of the chaining order is required since the opposite order is valid but does not lead to the same result.

CHAINING OF TRANSFORMATIONS BASED ON STATIC ANALYSIS

In this section, we present a static analysis based on formal definitions that allow us to determine whether a chain of transformations is valid. The definitions are given in Alloy, a lightweight specification language based on first order logic [START_REF] Jackson | Software Abstractions: Logic, Language, and Analysis[END_REF]. Our definitions are based on those in [START_REF] Mougenot | D-praxis: A peer-to-peer collaborative model editing framework[END_REF], which consider and encode a subset of MOF [START_REF] Omg | Meta Object Facility (MOF) 2.0 Query/View/Transformation Specification[END_REF].

Formalization of the concept of a model transformation

In [START_REF] Mougenot | D-praxis: A peer-to-peer collaborative model editing framework[END_REF], the authors consider models to be a finite set of model elements where each model element is typed by a metaclass and can own values for properties and references (complex concepts such as opposite references, derived property and subset values are not supported). In this section, we also use the syntactic domains defined in [START_REF] Mougenot | D-praxis: A peer-to-peer collaborative model editing framework[END_REF] and represented as sets of atoms: ME for model elements, MC for metaclasses, P for properties, R for references and V for property values. Listing 1 presents the Alloy signature of the considered domains. In Alloy, a set of atoms is specified by a signature.

Listing 1 alludes to the Alloy specification of models proposed in [START_REF] Mougenot | D-praxis: A peer-to-peer collaborative model editing framework[END_REF] Concerning the conformance relationship between model and metamodel, Listing 2 gives only the signature of the Alloy predicate. The body can be replaced by any definition provided in literature, such as the one specified in [START_REF] Paige | Metamodel-based model conformance and multiview consistency checking[END_REF].

1 pred conformsTo[m:Model, mm:MetaModel] {}
Listing 2: Conform Relationship.

Introduction of an intermediary representation

A transformation is defined in Listing 3 with its single input metamodel (mmin) and single output metamodel (mmout) (lines 2 and 3). To use alloy as a SAT solver, we consider that the execution of a transformation corresponds to the production of a single output model from one single input model (line 6). To that extent, we have specified within the transformation all of its input models (min) (line 4) and all of its output one (mout) (line 5); each input (resp. output) model conforming the input (resp. output) metamodel (line 9) (resp. (line 10)). Furthermore, we consider that the execution of a transformation might either have an impact (i.e., create, or modify, or delete), or read some elements of the input or output models. From these considerations, we introduce an intermediary representation of model transformation relying on the following properties: TImpacted, TCreated, TDeleted, TModified and TRead. Listing 4 introduces the corresponding domains that extend the abstract domain TProp. To each transformation is associated a set of properties (line 7 of Listing 3) specifying for which metaclasses, model elements can be impacted, (i.e., created, modified or deleted) or read. The metaclasses whose elements can be impacted are included in the set of the output metamodel metaclasses (line 11), whereas the metaclasses whose elements can be read, or impacted are included in the set of the input metamodel metaclasses (line 12). From these two lines it can be noticed that some metaclasses are present both in the input and the output metamodel.

Let The other properties are similarly built1 . They use predicates that define how elements are either deleted or modified. An instance of a metaclass mc has been deleted by the execution of the transformation if there exists, in the input model, some element whose metaclass is mc that belongs to the input model and not to the output model. A metaclass mc is considered modified if it exists at least one element of the input model whose mc is the metaclass that is changed i.e., its values or references in the input model are different from those in the output model.

Chaining Properties

Before specifying chaining constraints, let us review two definitions: chaining and commutativity (Listing 7). Two transformations t1 and t2 can be chained if the set of t1 output models is included in the set of t2 input models (lines 2 to 4). Two transformations t1 and t2 are commutative if t1 can be chained with t2, t2 can be chained with t1 and if each execution of the chain in one order or the other leads to the same result (lines 6 to 10). Thanks to Alloy and its SAT solver, we can state that two transformations t1 and t2 can be chained in either order, and lead to the same results on any input model, if the elements read or modified by one transformation are not impacted by the other, and if neither of the two transformations is the identity. Indeed, Alloy cannot find any example where the CommutativeTransformation predicate is true i.e., when some elements of the metamodels involved in t1 and in t2 are read or modified by one of the transformation and impacted by the other (line 7 of Listing 8); and none of the two transformations is the identity; and t1 and t2 are commutative. The static analysis based on Alloy leads to the following property (its symmetric is obviously true, t1 and t2 could be intertwined):

Property 1 (T Readt1 ∪ T M odif iedt1) ∩ T Impactedt2 = ∅ ⇔ t1 and t2 are not commutable.

Consider two transformations t1 and t2 where one or more metaclasses of the input metamodel also belong to the output metamodel of the transformation. If the intersection of the concepts read or modified by one transformation and the concepts impacted (i.e., whose some instances are created, deleted or modified) by the other is not empty, the two transformations are not commutable. In the other case, they are considered commutable.

IMPLEMENTATION

Alloy and its SAT solver enables us to formally define the commutative property. However, transformations are typically written in task-specific transformation languages such as QVTo or ATL, not Alloy, and so adaptations have to be performed. There are at least two alternatives: either encode the transformations in Alloy or adapt the formal static analysis to existing transformation languages. We chose the second alternative in order to preserve the transformations in the language used by the transformation developers. For this purpose, we have defined a metamodel corresponding to the intermediary representation and relying on the properties defined on Listings 5 and 6 in Alloy. The example transformations that we make use of in this paper are written in QVTo. In this section, we illustrate our static analysis using a transformation chain written in QVTo, to show the relationship between the intermediary representation and the languages that chain designers would apply.

Principle

Figure 1 sketches the approach. The intermediary representation is produced from a model representation 2 of the transformation to analyse. This step is itself implemented by a model transformation. Because it effectively takes a transformation as input, this is called a higher-order transformation (HOT) [START_REF] Bézivin | Model transformations? transformation models![END_REF]. Obviously, the choice of the language used to define this HOT is independent of the languages used to write the transformations that are to be chained. The static analysis we have presented is applied to models in the automatically generated intermediary representation. Following the analysis returns, a flag is raised to the chain designer. The chain designer can then verify if a chosen order is the one desired.

Creating the intermediary representation from a transformation

In the following paragraphs, we detail the way the TCreated, and TRead elements are produced from a QVTo transformation following the Alloy specifications of Section 3. The TDeleted and TModified are built similarly.

Creating elements

Three different ways are commonly used to create elements in the output model through a QVTo transformation: the mapping operation, the constructor operation and the object keyword. Each of these operations create a new instance (except inout mapping operation). This clearly respects the definition of TCreated and of the predicate ModelElement-InstanceOfMcHasBeenCreated specified in Section 3. Thus, for each mapping, constructor operation or object keyword contained in the transformation, a TCreated element referring the output metaclass is created.

Reading elements

In QVTo, read an instance is to access to it. Therefore, each time an existing instance is accessed and not modified (in order to conform to the TRead definition), a TRead referring to the involved metaclass is created. An instance can be accessed from different ways such as: from a mapping operation input, from successive navigations, from mapping operation body, etc.

Deleting and modifying elements

In QVTo, only the removeElement operation removes model elements instance of a given metaclass, from a model. Thus, for each removeElement operation, a TDeleted element, referring to the metaclass of the deleted element, is created. In QVTo, the inout mapping operation and the inout operation parameter may alter a property value or a reference of an element existing in the input and output models. Consequently, for each inout mapping or each inout parameter, a TModified element referring to the involved metaclass is created.

Set-based analysis

The set-based analysis checks if the elements read or modified by a transformation may be impacted (i.e. created, modified or deleted) by a second transformation. It corresponds to the implementation of the Alloy NonEmptyIntersect predicate and is concretely implemented through a QVTo transformation composed of a single query that checks the inclusion and returns a Boolean. This transformation also manages type inheritance. If in a metamodel, a metaclass X inherits from another metaclass Y and if an instance of X is created/modified/deleted by the transformation, then it also implicitly means than an instance of Y has been created/modified/deleted.

Similarly to what was done in Alloy in the previous section, this query is called twice. If one of the two intersections is not empty, the return value is true meaning that the two transformations are not commutative and thus conflicting. A flag is raised; the chain designer must verify that the chosen order is the one expected since we are in the case where both orders are valid (result from the type based analysis as described in [START_REF] Etien | Combining independent model transformations[END_REF]). In the opposite, if the returned value is false for both query calls, then whatever the models, the execution of t1 then t2 or t2 then t1 will lead to the same results. Until now the returned value (true or false) is printed in a log and is manually handled.

EVALUATION

In order to evaluate the efficiency of the proposed analysis (which we call TAction), we compared it with the inputoutput metamodels based analysis presented in [START_REF] Etien | Combining independent model transformations[END_REF].

For the experimentation, we considered three transformations sets. The first set is composed of 6 toy transformations, which correspond to elementary test cases. The second set is composed of 22 real transformations which are used to generate information systems. The 30 transformations of the third set also aim to generate information systems. The latter two sets have been provided by two different partner organisations, whereas the toy transformations were produced by the authors. For each set, we computed the number of commutable transformation pairs using the input-output metamodels based analysis of [START_REF] Etien | Combining independent model transformations[END_REF] 1 summarises, for each set, the number of transformations it contains; the number of transformation pairs and the number of commutable pairs identified for each of the two anslysis. For example, the second set can contain up to 253 commutable pairs. The input-output metamodels based analysis identified 33 pairs whereas our analysis only 26. For each set, the difference between the two analyses is very small (ranging from seven to five pairs).

We developed the TAction-based analysis to reduce false positives (i.e. pairs erroneously identified as commutative by the input-output metamodels based analysis). A deeper study of the results presented in Table 1 is required. The pairs considered commutable by the type based analysis were compared to the ones highlighted by the analysis proposed in this paper. Three cases occur: some pairs are identified by both analysis; some pairs only by the input-output metamodels based analysis (false positive), some pairs by only the TAction-based analysis (false negative). Some pairs are considered non commutative by the input-output metamodels based analysis, whereas the TAction-based analysis marked them as commutative. They result from the fact that the input-output metamodels based analysis takes into account the whole input and output metamodels associated to a transformation even if only a subset of concepts is used.

The false negatives disappear when the metamodels only contain metaclasses involved in the transformation (i.e. read or impacted). Furthermore, because of the over approximation due to the TAction-based analysis (e.g. dead code), and the used level of detail (i.e. based on metaclasses rather than properties) doubts persist concerning the results. In our examples, manual checking allows us to ensure that the obtained results are correct. . Furthermore, some pairs in these two cases were not considered as commutative whereas they are. The difference observed between set 2 and set 3, concerning false negative, arises from the adequacy of the metamodels to the transformations: the input-output metamodels based analysis uses metamodels with significant irrelevant concepts for the transformation. This evaluation highlights the addedvalue of our analysis based on a TAction-based analysis to check the commutativity between transformations: we can produce fewer false negatives.

Set

RELATED WORK

Several authors promote the decomposition of transformation into smaller ones in order to enhance their reusability and their maintainability. The way the chains are built from these small transformations varies according to the proposed approaches. Different operators exist: chaining [START_REF] Oldevik | Transformation composition modelling framework[END_REF][START_REF] Olsen | Aspects of reusable model transformations[END_REF], composition [START_REF] Cuadrado | Approaches for model transformation reuse: Factorization and composition[END_REF], that can eventually be conditional or parallel [START_REF] Rivera | Orchestrating ATL model transformations[END_REF], and loop [START_REF] Rivera | Orchestrating ATL model transformations[END_REF]. Mostly such operators can be applied when transformations are traditional heterogeneous ones and the chaining conditions relative to the inclusion of the output and input metamodels also [START_REF] Oldevik | Transformation composition modelling framework[END_REF][START_REF] Olsen | Aspects of reusable model transformations[END_REF][START_REF] Rivera | Orchestrating ATL model transformations[END_REF]. In [START_REF] Cuadrado | Approaches for model transformation reuse: Factorization and composition[END_REF], the transformations can be composed to create a new one but only if they have the same input metamodel; they may lead to a completely different output metamodel. In some approaches, for each small transformation, the required and provided concepts are explicitly manually identified by the chain designer using a profile [START_REF] Vanhooff | Breaking up the transformation chain[END_REF] or automatically identified based on the distinction between concepts copied and those mutated [START_REF] Chenouard | Automatically discovering hidden transformation chaining constraints[END_REF] or more generally using critical pair analysis [START_REF] Lambers | Sufficient criteria for applicability and non-applicability of rule sequences[END_REF][START_REF] Ehrig | Fundamentals of Algebraic Graph Transformation[END_REF][START_REF] Mens | Detecting Structural Refactoring Conflicts Using Critical Pair Analysis[END_REF]. However, in these latter cases, they only deal with endogenous transformations (even if they suggest that an extension to heterogeneous transformations is possible) and do not care about commutativity.

In [START_REF] Lambers | Sufficient criteria for applicability and non-applicability of rule sequences[END_REF] and [START_REF] Ehrig | Fundamentals of Algebraic Graph Transformation[END_REF], the authors also use graph theory to precisely specify when two transformation rules of pure refactoring are conflicting. If the critical pair analysis is adapted to the detection of conflicts between rules, it cannot be adapted between transformations. Critical pair analysis is based on the single pushout semantics where rules have a left-hand side, a right-hand side, and eventually a NAC. Such a semantics is not general enough to considered a whole transformation expressed in different transformation languages. Our relies on the same intuition than the one underlying the parallel and sequential independence described in [START_REF] Ehrig | Fundamentals of Algebraic Graph Transformation[END_REF]. Nevertheless it takes into account the transformation in its whole and not each of its rules. To that extent, it proposes the introduction of an intermediary representation to detect conflict between transformations.

In this paper we provide a static analysis based on an intermediary representation where concepts created, modified or read are distinguished. The designer has nothing to specify by hand, everything is automatically computed. Our approach does not impose any constraint on the metamodels of the transformation. It takes all its benefits when the input and output of the transformations are different but contain some metaclasses in common.

CONCLUSION

In this paper, we have proposed an approach to verify the chaining of model transformations. Our mechanism aims at checking if two transformations are commutative, i.e. if they can be chained in either order and deliver the same results no matter the input. This mechanism relaxes the constraints previously established by distinguishing the read and impacted (i.e. created, deleted or modified) concepts.

We formally proposed a static analysis relying on an intermediary representation independent from any transformation language. We have implemented our approach on top of QVTo in order to check the commutativity of transformations written with this language. If the two studied transformations are not commutative, a flag is raised to the chain developer in order she verifies if the chosen order is the expected one. As any other static analysis, our approach suffers from over approximations due to dead code for example. We also foresee to take into account not only the metaclasses but also the properties. With our approach two transformations modifying two different properties of the same metaclass are not currently considered as commutative whereas they can be.

As a conclusion, our investigations have been validated using QVTo but do not depend on any model transformation technique. The model transformations are handled as functions having one input and one output metamodel possibly sharing some concepts. This approach refines the existing chaining constraints and will help to design chains from existing transformations.

16 } 1 :

 161 Listing Signature of Model and Metamodel.

11 ( 13 }

 1113 min | conformsTo[mi,mmin] 10 all mo:mout | conformsTo[mo,mmout] prop).TImpacted in mmout.mc 12 (prop).TRead + TImpacted in mmin.mc Listing 3: Transformation.

1 abstract sig TProp {} 2 sig

 12 TImpacted, TCreated, TModified, TDeleted, TRead extends TProp {} Listing 4: Properties associated to transformation.

14 } 15 } 6 :

 14156 Listing Definition of the used predicates.

1 //Chaining 2 pred 4 } 5 6//Commutativity 7 pred 11 } 7 :

 12457117 ChainingTransformation[t1,t2:Transformation] { 3 t1.mout in t2.min Commutativy[t1,t2:Transformation] { 8 ChainingTransformation[t1,t2] 9 ChainingTransformation[t2,t1] 10 all m:t1.min+t2.min | t2.execution.(t1.execution.m) = t1 .execution.(t2.execution.m)Listing Chaining and Commutativity properties.

1 pred 14 } 8 :

 1148 NonEmptyIntersect[t1,t2:Transformation] { 2 some ((t1.prop).TRead + (t1.prop).TModified )& ((t2.prop ).TImpacted)Listing Commutative Transformations checking properties.

Figure 1 :

 1 Figure 1: Sketch of the approach 2 Most transformation tools such as QVTo or ATL provide such a representation.

  us now precisely specify to what these properties correspond. A metaclass belongs to the TCreated set, if there exists some execution of the transformation that creates a model element instance of this metaclass. The TCreated property is specified by the TransformationCreated-Compliant predicate (Listing 5 lines 11 to 19) that uses the predicates ModelElementInstanceOfMcHasBeenCreated, which defines the meaning of instance creation.

	1	//An instance has been created (Used to define TCreated)
	2	pred ModelElementInstanceOfMcHasBeenCreated[mi:Model , mo
		:Model , mc:MC] {
	6	}
	7	}
	8
	9	//TCreated
	8	}
	9
	10	//All instances are unchanged (Used to define TRead)
	11	pred ModelElementAreReadOnly[mi:Model, mo:Model , mc:MC]
		{
	12	all melt:mi.me {
	13	melt.(mi.class) = mc implies ModelElementIsUnchanged[mi,
		mo,melt]

3 some melt:mo.me { 4 melt.(mo.class) = mc 5 melt !in mi.me 10 pred TransformationCreatedCompliant[t:Transformation] { 11 all cmc: (t.prop).TCreated { 12 some mi,mo: Model { 13 (mi->mo) in t.execution 14 ModelElementInstanceOfMcHasBeenCreated[mi,mo, cmc] 15 } 16 } 17 } 18 19 pred TransformationCreatedConstrained[t:Transformation] { 20 some (t.prop).TCreated 21 } Listing 5: Details of the TCreated property. A metaclass mc is considered read only if each element of the input model whose mc is the metaclass, is unchanged (Listing 6 lines 11 to 15). A model element is unchanged if it belongs to the input and the output models and if its metaclass, its property values and its reference values are the same in both models (lines 2 to 8). 1 //The model element remains unchanged (Used to define TRead) 2 pred ModelElementIsUnchanged[mi:Model, mo:Model,melt:ME] { 3 melt in mi.me 4 melt in mo.me 5 melt.(mi.class) = melt.(mo.class) 6 melt.(mi.valueP) = melt.(mo.valueP) 7 melt.(mi.valueR) = melt.(mo.valueR)

Table 1 :

 1 and our TAction-based analysis. Commutativity results

	Set #Transfo. Pair in.out. MM	TActions
	1	6	21	13 (61.90%)	6 (38.06%)
	2	22	253	33 (13.04%)	26 (10.27%)
	3	30	465	45 (9.68%)	40 (8.60%)
	Table

Table 2 :

 2 Numbers of errors and approximations using TAction-based analysisTable2gathers the figures for each set computed by the TAction-based analysis. Once again, set 1 has been developed for test purposes and the resulting figures are not trustworthy. In contrast, in the two other sets, around 30% of the commutative pairs identified by the input-output metamodels based analysis are in fact not commutative (false positive)

		Common	False positive False negative
	1	5 (38.46%)	8 (61.53%)	1 (7.69%)
	2	23 (69.69%)	10 (30.30%)	3 (9.09%)
	3	30 (66.66%)	15 (33.33%)	10 (22.22%)

The complete code can be found at http://www.lifl.fr/∼etien/commutativity


Jean-Paul Chehab

Denys Dutykh
email: denys.dutykh@univ-savoie.fr

J.-P Chehab

On time relaxed schemes and formulations for dispersive wave equations

Keywords: dispersive wave equations, shallow water flows, relaxation, quasi-compressibility. MSC: [2010] 74J15 (primary), 74S10, 74J30 (secondary) PACS: [2010] 47.35.Bb (primary), 47.35.Fg, 47.85.Dh (secondary) dispersive wave equations, shallow water flows, relaxation, quasicompressibility

The numerical simulation of nonlinear dispersive waves is a central research topic of many investigations in the nonlinear wave community. Simple and robust solvers are needed for numerical studies of water waves as well. The main difficulties arise in the numerical approximation of high order derivatives and in severe stability restrictions on the time step, when explicit schemes are used. In this study we propose new relaxed system formulations which approximate the initial dispersive wave equation. However, the resulting relaxed system involves first order derivatives only and it is written in the form of an evolution problem. Thus, many standard methods can be applied to solve the relaxed problem numerically. In this article we illustrate the application of the new relaxed scheme on the classical Korteweg-de Vries equation as a prototype of stiff dispersive PDEs.

Introduction

Nonlinear dispersive waves arise in various fields of science including the solid [START_REF] Porubov | Propagation of localized longitudinal strain waves in a plate in the presence of cubic nonlinearity[END_REF] and fluid mechanics [20]. Below we focus on models stemming from the free surface hydrodynamics mainly due to scientific interests of the authors. In this field, analytical solutions (such as solitary or cnoidal waves [8,10]) are available only for some simplified models. That is why the numerical simulation remains one of the main tools in nonlinear dispersive wave studies. In this work the authors have an ambition to propose a novel strategy to tackle numerically these problems.

Let us review first the modern numerical approaches to this problem. Historically, the finite differences were applied to dispersive wave problems [START_REF] Ismail | A finite difference method for Korteweg-de Vries like equation with nonlinear dispersion[END_REF]. Today, the trend is to apply Galerkin-type methods (e.g. FEM) in smooth situations [START_REF] Antonopoulos | Galerkin approximations of the periodic solutions of Boussinesq systems[END_REF][START_REF] Mitsotakis | On the Galerkin/Finite-Element Method for the Serre Equations[END_REF] and finite volumes [START_REF] Dutykh | Finite volume schemes for dispersive wave propagation and runup[END_REF]15] or coupled finite volumes (for the hyperbolic part)/finite differences (for the dispersive part) via the operator splitting [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF]. Notice that there is an effort to develop discontinuous Galerkin-type methods for dispersive wave equations as well [16,[START_REF] Yan | A local discontinuous Galerkin method for KdV type equations[END_REF], which combine the robustness of finite volumes with the accuracy of finite elements. However, the 'price to pay' (i.e. the CPU time) remains excessively high turning these methods into a luxury limousine rather than a working horse of numerical analysis. Finally, for periodic situations one can apply highly efficient Fourier-type pseudo-spectral methods [START_REF] Dutykh | Finite volume and pseudo-spectral schemes for the fully nonlinear 1D Serre equations[END_REF]. On the borderline between finite differences and pseudo-spectral methods there exist compact finite difference schemes proposed by [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] with spectral-like resolution. To give an example, these schemes were applied to the Serre equations in [9].

In the present study we propose to change the numerical strategy in contrast to studies described above. The main idea consists in modifying the governing equations by slightly perturbing them with some ad-hoc terms. In this way, we gain the structure suitable for numerical simulations and we hope that solutions of the perturbed and un-perturbed problems will remain close provided that perturbation parameter is small. This idea is inspired by pseudo-compressibility methods proposed for the numerical simulation of incompressible Navier-Stokes equations (see e.g. [START_REF] Kameyama | Multigrid iterative algorithm using pseudocompressibility for three-dimensional mantle convection with strongly variable viscosity[END_REF]). The main philosophy of this approach is as follows. When we solve numerically a differential equation, some errors are introduced due to the underlying discretization process. So, perhaps, for the sake of convenience, one can perturb also the governing equations within the same discretization error which is already introduced. And if everything is done judiciously, the end user will not even see the difference between the relaxed and the original problems solutions. It makes a lot of 'ifs', but nevertheless we undertake this programme below for some well-known dispersive wave equations. We would like to mention here also another relaxation technique based on the application of local time-averaging operators [2]. A relaxation scheme for the KdV equation was proposed in [START_REF] Benkhaldoun | New finite-volume relaxation methods for the third-order differential equations[END_REF]. However, their formulation contains 2 nd order derivatives. Below, we will propose an alternative formulation which has an advantage to involve only the 1 st order derivatives. The idea we employ can be traced back at least to various local * formulations used in continuous [START_REF] Walkley | A finite element method for the one-dimensional extended Boussinesq equations[END_REF][START_REF] Walkley | A finite element method for the two-dimensional extended Boussinesq equations[END_REF] and discontinuous [START_REF] Levy | Local discontinuous Galerkin methods for nonlinear dispersive equations[END_REF][START_REF] Yan | A local discontinuous Galerkin method for KdV type equations[END_REF] Galerkin methods. However, we push it one step further towards the so-called relaxed local formulations.

The present study is organized as follows. In Section 2 we propose the relaxed formulations for several well-known dispersive wave models. Some mathematical properties of a relaxed formulation for the KdV equation are discussed in Section 3. Our numerical approach based on this relaxed formulation of the KdV equation is outlined in Section 4. The first numerical results are presented in the same Section. Finally, in Section 5 we outline the main conclusions and perspectives of our study. The relaxation of some other dispersive wave equations is discussed in Appendix A.

Mathematical models

Instead of working out the most general situation, in the present study we follow the socalled Gelfand † principle which states that a theory should be illustrated on the simplest non-trivial example. Below we consider one such example and two others are treated in Appendix A. All cases which steem from the water wave theory [START_REF] Stoker | Water Waves: The mathematical theory with applications[END_REF]. For each model we develop the corresponding relaxed formulation.

Korteweg-de Vries equation

The celebrated Korteweg-de Vries (KdV) equation was derived independently by J. Boussinesq [7] and D. Korteweg & G. de Vries [START_REF] Korteweg | On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves[END_REF] at the end of the XIX th century and in the scaled form it reads:

u t + u u x + u xxx = 0 .
(2.1)

The variable u (x, t) may be interpreted physically as the free surface elevation or the horizontal depth-integrated fluid velocity. There exist transformations in order to obtain the canonical form (4.1) [20]. We shall rewrite the last equation in the conservative form in order to obtain the flux of the evolution variable u :

u t + 1 2 u 2 + u xx x = 0 . (2.2)
Now, the order of derivatives has to be lowered. It can be done similarly to the ODE theory by introducing additional variables:

u t + 1 2 u 2 + w x = 0 , u x = v , v x = w .
At this step we are compatible with local formulations used in the framework of discontinuous Galerkin methods [START_REF] Levy | Local discontinuous Galerkin methods for nonlinear dispersive equations[END_REF][START_REF] Yan | A local discontinuous Galerkin method for KdV type equations[END_REF]. It is obvious that the last system is completely equivalent to equation (2.2) (at least for the smooth solutions). However, the last system is not of the evolution type. It is here that the relaxation comes into the play in order to correct this shortcoming:

u t + 1 2 u 2 + w x = 0 , (2.3) δ v t + v -u x = 0 , (2.4
)

δ w t + w -v x = 0 , (2.5)
where δ ≪ 1 is a small parameter. It can be seen that in the limit δ → 0 we recover the original formulation (4.1) (at least formally). Equations (2.4), (2.5) allow also to assess the physical meaning of the parameter δ . Indeed, the terms v and δ v t (and w with δ w t ) need to have the same units. Thus, the multiplication by δ has to compensate the derivative with respect to time. Henceforth, from dimensional arguments, δ has the unit of time.

Properties

Above we presented relaxed formulations for three well-known equations. However, the properties of new formulations have to be studied. From now on we focus on the KdV equation only (2.3) -(2.5) in order to understand the relaxation effects exhaustively. Other equations will be tackled in following studies.

Dispersion relation analysis

In order to study the dispersive properties, first we have to linearize equations (2.3) -(2.5):

u t + w x = 0 , (3.1)
δ v t + v -u x = 0 , (3.2) δ w t + w -v x = 0 . (3.3)
Now, we look for plane-wave solutions of the form

u (x, t) = u 0 • e i (k x -ω t) , v (x, t) = v 0 • e i (k x -ω t) , w (x, t) = w 0 • e i (k x -ω t) ,
where k is the wavenumber, ω is the frequency and u 0 , v 0 , w 0 ∈ R are some real amplitudes. By substituting * the plane-wave ansatzä into a linear system (3.1) -(3.3), we obtain a system of linear algebraic equations, where unknowns are u 0 , v 0 , w 0 :

-i ω u 0 + i k w 0 = 0 , -i δ ω v 0 + v 0 -i k u 0 = 0 , -i δ ω w 0 + w 0 -i k v 0 = 0 .
The necessary condition to have non-trivial solutions to the above linear system of equations is that its determinant is equal to zero:

det -i ω 0 i k -i k 1 -i δ ω 0 0 -i k 1 -i δ ω = 0 .
After some simplifications we obtain the desired dispersion relation between ω and k :

ω + k 3 -2 i δ ω 2 -δ 2 ω 3 = 0 . (3.4)
So, one can see that by taking the limit δ → 0 in the last equation, we recover the classical KdV dispersion relation ω = -k 3 . This polynomial equation (3.4) in wave frequency ω admits in general three complex roots. Their general expression is rather cumbersome. However, we can construct their asymptotic expansions:

ω 1 (k) = -k 3 + 2 i δ k 6 + 7 k 9 δ 2 + O (δ 3 ) , (3.5) ω 2, 3 (k) = - i δ ± i √ -i k 3/2 √ δ + 1 2 k 3 ± 5 8 √ -i k 9/2 √ δ -i k 6 δ + O (δ 3/2 ) . (3.6)
These expansions (3.5), (3.6) might be used to interpret the behaviour of numerical solutions. Since we deal with singular perturbations in equation (3.4), the passage to the limit δ → 0 in solutions ω 2, 3 (k) is far from being trivial. However, from the expansion (3.5) we learn that parameter δ plays the both rôles of dissipation and dispersion since it enters into real and imaginary parts of the regular branch ω 1 (k) . We would like to highlight the fact that the relaxation parameter δ affects the dispersive properties of the regular branch (i.e. Re ω 1 (k)) only to the second order in δ . Thus, we can conclude that the proposed relaxation approach is very careful with respect to the dispersion relation of the original KdV equation. Remark 1. Taking the limit in the polynomial equation (3.4) turns out to be an easy task, as it often happens, than taking the limit δ → 0 in the solutions ω 2, 3 (k) . On the other hand, the branch ω 1 (k) seems to be completely regular from this point of view. Two other branches ω 2, 3 (k) are artificial modes introduced by relaxation.

Numerical methods and results

In this Section we are going to study numerically the relaxed formulations. For this purpose we propose three different discretisations of the KdV equation. One of them is based on the relaxed formulation proposed above and two others are popular finite difference schemes (Crank-Nicolson and Sanz-Serna) based on the original KdV equation. The numerical advantages of the relaxed formulation are illustrated below.

Before considering the time schemes, we focus on the discretization in space which is realized with finite differences. The use of compact schemes as presented in [23] allows to reach a high order of accuracy, with a spectral-like resolution, while preserving the capabilities of the finite differences approach, including the implementation on a general Cartesian grids and for a variety of boundary conditions. When dealing with surface wave equations, the high level of accuracy is an important property to capture correctly the active frequencies of the solution without working on very fine grids. Hence, the computation can be done with a reasonable computational time.

At first we recall briefly the principle of the Compact Schemes (CS) and we give the discretization matrices that will be used for the simulation of the models presented in Section 2.

In two words, the compact schemes consist in approaching a linear operator (differentiation as well as interpolation) by a rational (instead of polynomial-like) finite differences scheme. These schemes are then implicit, this allows to increase the accuracy and mimic the spectral global dependence, see [23] for more details. The approximations at grid points of the operators applied to a regular function u is realized as follows.

Let U def := (U 1 , . . . , U N ) ⊤ be a vector whose the components are the approximations u at (regularly spaced) grid points

x i = i h , i = 1, . . . , N , here h = 1 N
is the spatial step-size. We compute approximations of V i = L (u) (x i ) as solution of a system

P • V = Q U ,
the approximation matrix is then formally B = P -1 • Q . For simplicity we present here fourth order accurate CS. When considering periodic boundary conditions, the matrices P and Q are for the first order derivative in space:

P def :=       1 1 4 1 4 1 4 1 . . . . . . . . . 1 4 1 4 1 4 1       , Q def := 1 2 h         0 3 2 -3 2 -3 2 0 3 2 . . . . . . . . . -3 2 0 3 2 3 2 -3 2 0         ,
For the second order derivative in space we have the following pair of matrices:

P =        
1         , Q = 1 h 2             12
            .
Finally, for the third order derivative in space we have

P =         1 1 2 1 2 1 2 1 1 2 . . . . . . . . . 1 2 1 1 2 1 2 1 2 1         , Q = 2 h 3         0 -1 1 2 -1 2 1 1 0 -1 1 2 -1 2 . . . . . . . . . 1 2 -1 2 1 0 -1 -1 1 2 -1 2 1 0        
.

We now present the semi-implicit relaxed time scheme we will use for the simulation, stressing out its enhanced stability as compared to classical semi-implicit schemes (Crank-Nicolson's for the linear terms and forward Euler's for the nonlinear ones) while giving comparable results to the ones obtained by fully nonlinear Sanz-Serna's which is second order accurate in time and unconditionally stable.

Let D x , D xx and D xxx be the discretization matrices of the first, second and third order derivative on a N regularly spaced grid points, with periodic boundary conditions; these matrices will be obtained in practice with compact schemes as presented above; I N will denote the N × N identity matrix. The relaxed scheme produces the iterations

U (k+1) -U (k) ∆t + 1 2 D x (W (k+1) + W (k) ) + (U (k) ) 2 = 0 , δ V (k+1) -V (k) ∆t + 1 2 (V (k+1) + V (k) ) -D x (U (k+1) + U (k) ) = 0 , δ W (k+1) -W (k) ∆t + 1 2 (W (k+1) + W (k) ) -D x (V (k+1) + V (k) ) = 0 .
We rewrite this coupled system by blocs as

   I N 0 ∆t 2 -∆t 2 D x (δ + ∆t 2 ) I N 0 0 -∆t 2 D x (δ + ∆t 2 ) I N    •    U (k+1 ) V (k+1 ) W (k+1 )    =    I N 0 -∆t 2 ∆t 2 D x (δ -∆t 2 ) I N 0 0 ∆t 2 D x (δ -∆t 2 ) I N    •    U (k) V (k) W (k)    +    -∆t 2 D x (U (k) ) 2 0 0    .
We can now introduce the iteration matrix

M r def := (M i r ) -1 • M e r =    I N 0 ∆t 2 -∆t 2 D x (δ + ∆t 2 ) I N 0 0 -∆t 2 D x (δ + ∆t 2 ) I N    -1 •    I N 0 -∆t 2 ∆t 2 D x (δ -∆t 2 ) I N 0 0 ∆t 2 D x (δ -∆t 2 ) I N    .
We can resume the relaxed scheme as follows Algorithm 1 Relaxed Scheme for the KdV equation.

1: U (0) is given 2: Set V (0) = D x U (0) , W (0) = D x V (0) 3: Set Z (0) = U (0) , V (0) , W (0) ⊤ 4: for k = 0, 1, . . . do 5: Set F (k) = -∆t 2 D x (U (k) ) 2 , 0, 0 ⊤ 6: Solve M i r • Z (k + 1) = M e r • Z (k) + F (k) 7: Set U (k+1) = Z (k + 1) (1 : N) , V (k + 1) = Z (k + 1) (N + 1 : 2 N) , W (k + 1) = Z (k + 1) (2 N + 1 : 3 N) 8: end for
We now present the two references schemes to which we will compare the relaxed one. We set

M i CN def := I N + ∆t 2 D x x x , M e CN def := I N - ∆t 2 D x x x
and we define the corresponding iteration matrix M CN as

M CN def := (M i CN ) -1 • M e CN .
Algorithm 2 Fully nonlinear Sanz-Serna's scheme.

1: U (0) is given 2: 3: for k = 0, 1, . . . do 4: Solve M i CN U (k + 1) = M e CN U (k) -∆t 2 D x U (k + 1) + U (k) 2 2 5: end for
We will use also

Linear stability analysis

Before comparing the stability properties of the classical CN and of the relaxed schemes, we give hereafter a simple but instructive result illustrating the advantage of the compact schemes in a context in which the spectral properties of the operators must be restituted Algorithm 3 Semi-implicit Crank-Nicolson scheme.

1: U (0) is given 2: 3: for k = 0, 1, . . . do 4: Solve M i CN U (k + 1) = M e CN U (k) -∆t 2 D x U (k)
û (k + 1) m - û (k) m ∆t - i 2 m 3 û (k + 1) m + û (k) m = 0 , hence û (k + 1) m = 1 + i ∆t 2 m 3 1 -i ∆t 2 m 3 û (k) m .
The Fourier symbol of the CN operator is

1 + i ∆t 2 m 3 1 -i ∆t 2 m 3
and its values are all displayed on the unit circle, for any ∆t > 0 . It is then an important property to be captured by the spectrum of the Crank-Nicolson matrix, when dealing with finite differences. We give hereafter in Figure 1 the spectra of the CN matrix when the 2 nd order FD and the 4 th order compact FD are used. We observe that the compact scheme allows to capture numerically the spectral properties of the linear operator as in the analytical Fourier analysis, while the 2 nd order finite differences succeed to capture correctly only a small part of the spectrum. Comparable results are obtained in Figure 2 when considering 6 th order compact FD: they illustrate the capabilities of the compact schemes to mimmic the spectral properties of the linear propagation operator.

Enhanced stability analysis

A first way to understand the effect of the relaxation in the stability of the time marching schemes is to consider the linear system (3.1). We have used the 6 th order compact schemes for the discretization of the spatial derivatives, as presented in the beginning of Section 4. We compare here in the complex plane the eigenvalues of the two matrices used in the numerical schemes: M r and M CN . We recall that these matrices are not of the same size so there is no reason that their relative spectra coincide. The eigenvalues of M CN are mapped on the unit circle (see Figure 1), we present hereafter comparison with those of M r for different values of N , ∆t and δ (the relaxation parameter). At first, we take a very small value of δ = 10 -17 , however the results we obtain are identical all for larger but still small values of δ , up to 10 -5 , as illustrated in Figure 3. For larger values of δ , say δ = 10 -4 , δ = 10 -3 , δ = 10 -1 , eigenvalues of M r are placed outside the unit disk making the relaxed scheme unstable, see Figures 4, 5 and6. Below, in Table 1, we give the maximum and the minimum of the modulus of the eigenvalues of the M r and M CN matrices for different values of the time step ∆t and relaxation parameter δ . We can observe the influence of δ on the spectral radius of M r : the relaxed scheme becomes unstable, i.e. ρ (M r ) > 1 , for not sufficiently small values of δ , say, e.g. δ 10 -4 ; for small values of δ , we have ρ (M r ) = 1 , the eigenvalues of M r are perfectly matched on the unit circle and the scheme is then unconditionally stable.

Simulation of a soliton

To illustrate the effect of the time relaxation, we consider the simulation of a soliton on the domain I = 0, ℓ , for different values of δ and of ∆t , namely

u (x, t) = a sech 2 1 2 κ (x -x 0 -c t) , (4.1)
with c = a 3 and κ = √ c . In our experiments presented below we use a = 0.8 , ℓ = 100 and x 0 = ℓ 2 . The interval I is discretized with N = 400 equidistant points and the space discretization is realized by using 6 th compact schemes. Of course, due to the nonlinear orbital stability properties of the solitons [START_REF] Weinstein | Existence and dynamic stability of solitary wave solutions of equations arising in long wave propagation[END_REF], the approximation to the solitary wave by the numerical solution of the relaxed system cannot be considered in a long time interval, making the validation delicate. Indeed, the classical error measures such as L p norms are irrelevant for orbits, meanwhile a numerical method approximates a classical solution and not an orbit. We observe that the relaxation allows to approach the exact solution with an expected level of accuracy. For instance, when δ = 10 -17 , the solutions computed by the different schemes coincide on a fairly long time interval. With δ = 10 -4 and ∆t = 0.1 , the solution coincide on a smaller time interval as it is expected. More generally, we notice in our numerical experiments that the pairs (δ , ∆t) that make the relaxed scheme stable in the linear case are operant for the nonlinear KdV equation as well * . The contrary is also * However, we do not have a mathematical proof for this statement. observed, a numerical instability holds, for example, for ∆t = 0.005 and δ = 0.001 , see Table 1. The numerical result for the soliton propagation are shown in Figure 7 for two values of the relaxation parameter δ = 10 -17 and δ = 10 -4 . On the right panels we show the evolution of the L 2 error computed thanks to the exact Solution (4.1).

Second order relaxed times schemes

As illustrates above, relaxed schemes allow a correct numerical integration using only approximation to ∂ x which is generally easily available in a number of situations, while the approximation to ∂ x x x can be very tricky to build. However, the relaxation scheme is only first order accurate in time and as the semi-implicit Euler scheme, it reveals not to be suited for long time interval simulations. To overcome this drawback we propose here to reach a second order of accuracy by using a classical Richardson extrapolation. In two words, the numerical time integration of ODE

du dt = F (u) ,
by the forward Euler scheme defines the iterations Eigenvalues Iter. Mat.

u k + 1 = u k + ∆t F (u k ) = G ∆t (u k ) ,
relaxed scheme CN scheme

(b) ∆t = 10 -2 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3
Eigenvalues Iter. Mat. which are first order accurate approximations to u (k ∆t) . The Richardson extrapolated sequence is defined by

v 1 = G ∆t (u k ) , v 2 , 0 = G ∆t/2 (u k ) , v 2 , 1 = G ∆t/2 (v 2 , 0 ) , u k + 1 = 2 v 2 , 1 -v 1 ,
and is second order accurate in time. We will start here from a simple IMEX method, says backward Euler for the linear terms and forward for the nonlinear ones: in other words if F (u) writes as F (u) = -A u + H (u) , then the propagator is formally defined as

G ∆t (u k ) = (I + ∆t A) -1 • u k + ∆t H (u k ) .
Table 1. Minimal and maximal values of eigenvalues of iteration matrices for N = 400 and ℓ = 100 . Here, 6 th order Compact Schemes are used.

N ∆t δ M CN : σ CN , ρ CN M r : σ r , ρ r 10 -3 10 -17 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1 10 -2 10 -17 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1 10 -1 10 -17 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1 5 × 10 -1 10 -17 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1 10 -3 10 -10 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.000000047 10 -2 10 -10 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.00000007 10 -1 10 -10 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.000000008 5 × 10 -1 10 -10 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.0000000015 10 -3 10 -8 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.000004775 10 -2 10 -8 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.00000691 10 -1 10 -8 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.000000798 5 × 10 -1 10 -8 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.0000001599 10 -3 10 -4 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.047217 10 -2 10 -4 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.070529 10 -1 10 -4 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.00796 5 × 10 -1 10 -4 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.001598 10 -3 10 -3 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.178956 10 -2 10 -3 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.56593 10 -1 10 -3 σ CN = 1, ρ CN = 1 σ r = 1 , ρ r = 1.07633 5 × 10 -1 10 -3 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.015844 10 -3 10 -1 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.02591 10 -2 10 -1 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.290614 10 -1 10 -1 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 2.90107786 5 × 10 -1 10 -1 σ CN = 1, ρ CN = 1 σ r = 1, ρ r = 1.909388
The IMEX relaxed scheme consists then in solving the following linear system at each step: The matrix of the system is noted by M (IMEX) ∆t , δ

   I N 0 -∆t D x -∆t D x (δ + ∆t) I N 0 0 -∆t D x (δ + ∆t) I N    •    U (k + 1 ) V (k + 1 ) W (k + 1 )    =    U (k) V (k) W (k)    +    -∆t 2 D x (U (k) ) 2 0 0    .
. The Extrapolated Relaxed schemes writes as Algorithm 4 : Extrapolated Relaxed Scheme

1: u (0) given 2: 3: for k = 0, 1, • • • until convergence do 4: Solve M (IMEX) ∆t/2, δ v 1 = -∆t 2 F (u (k ) , 5: Solve M (IMEX) ∆t/2, δ v 2 = -∆t 2 F (u 1 ) , 6: Solve M (IMEX) ∆t, δ v 3 = -∆t F (u (k) ) , 7:
Set u (k + 1) = 2 u 2 -u 3 . 8: end for Hereafter, in Figure 8 -9, we present the comparison of the evolution of the soliton and its numerical approximations (Sanz-Serna and extrapolated Relaxed and IMEX schemes). We observe that the time extrapolation allows to approach the exact solution with a good accuracy on longer time intervals, when considering different values of δ and ∆t ; this is notable particularly, e.g. when taking δ = 10 -4 and ∆t = 10 -1 .

Discussion

Above we presented some rationale behind time-relaxed formulations for several wellknown dispersive wave equations. The main conclusions and perspectives of our study are outlined below.

Conclusions

In this study we presented an approximate reformulation for several dispersive wave equations. This formulation was inspired somehow by quasi-(or pseudo-) compressibility methods to solve incompressible Navier-Stokes equations [START_REF] Kameyama | Multigrid iterative algorithm using pseudocompressibility for three-dimensional mantle convection with strongly variable viscosity[END_REF]. So, by following this 'philosophy' we proposed relaxed formulations for celebrated Korteweg-de Vries (KdV), Benjamin-Bona-Mahony (BBM) and Peregrine's system of equations. However, it is obvious that the same technique can be extended to many other scalar and vectorial dispersive wave equations. This formulation has a simple advantage to involve first order derivatives only while being in the form of a coupled evolution problem. This is the main difference with various local (or modified) reformulations used in continuous and discontinuous Galerkin methods [START_REF] Levy | Local discontinuous Galerkin methods for nonlinear dispersive equations[END_REF][START_REF] Walkley | A finite element method for the one-dimensional extended Boussinesq equations[END_REF][START_REF] Walkley | A finite element method for the two-dimensional extended Boussinesq equations[END_REF][START_REF] Yan | A local discontinuous Galerkin method for KdV type equations[END_REF]. Many standard numerical methods can be applied to solve the proposed relaxed formulation numerically. In the present study we applied compact finite differences (with spectral-like resolution) to discretize the problem in space along with a simple time stepping. The presented numerical tests and validations for the KdV equation show that the relaxed discrete formulation possesses a larger CFL-type stability limit comparing to similar compact discretizations of the standard KdV equation (without relaxation). This preliminary conclusion indicates that relaxed schemes might be good candidates for the numerical simulation of notoriously stiff systems of equations such as the KdV-KdV Boussinesq-type system extensively studied numerically in e.g. [5].

Perspectives

Above we presented some numerical illustrations for the classical KdV equation only. As it was mentioned in the previous Section, this formulation was extended to other weaklynonlinear models as well. The main goal consists in extending these time-relaxation numerical methods to Boussinesq-type systems of weakly nonlinear weakly dispersive equations such as the celebrated classical Peregrine [START_REF] Peregrine | Long waves on a beach[END_REF], modified Peregrine [START_REF] Durán | Peregrine's System Revisited[END_REF] or even fully nonlinear Serre-Green-Naghdi (SGN) equations [START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF] (some more conventional numerical strategies for these equations were outlined in e.g. [START_REF] Dutykh | Finite volume and pseudo-spectral schemes for the fully nonlinear 1D Serre equations[END_REF][START_REF] Mitsotakis | On the nonlinear dynamics of the traveling-wave solutions of the Serre system[END_REF][START_REF] Mitsotakis | On the Galerkin/Finite-Element Method for the Serre Equations[END_REF]).

On a more theoretical side, we would like to obtain the estimations of the difference between the perturbed and unperturbed problems solutions. In other words, we would like to have theoretical arguments to state that the proposed relaxation process generates solutions which remain close to those of the original equation. This is the main point on our 'theoretical' agenda.

2 5 :

 5 end for at the discrete level. Writing the CN scheme for the Airy equation, we have the following induction relation among Fourier coefficients attached to the m th frequency

FourierFigure 1 .

 1 Figure1. Spectra of the M CN matrices for the 2 nd (red markers) and 4 th order (black markers) finite differences schemes: ℓ = 20 , N = 100 , ∆t = 10 -3 and ∆t = 10 -2 .

-

Figure 2 .

 2 Figure2. Spectra of the M CN matrices for the 2 nd (red markers) and 6 th order (black markers) finite differences schemes: ℓ = 100 , N = 400 , ∆t = 10 -3 and ∆t = 10 -2 .

1 Figure 3 .

 13 Figure 3. Spectra of M r and M CN matrices. For the spatial discretisation we use ℓ = 100, N = 400 points and the relaxation parameter δ = 10 -17 . Various values of the time step are taken.

  ∆t = 5 × 10 -1

Figure 4 .

 4 Figure 4. Spectra of M r and M CN matrices. For the spatial discretisation we use ℓ = 100, N = 400 points and the relaxation parameter δ = 10 -4 . Various values of the time step are taken.

1 Figure 5 .

 15 Figure 5. Spectra of M r and M CN matrices. For the spatial discretisation we use ℓ = 100, N = 400 points and the relaxation parameter δ = 10 -3 . Various values of the time step are taken.

1 Figure 6 .

 16 Figure 6. Spectra of M r and M CN matrices. For the spatial discretisation we use ℓ = 100, N = 400 points and the relaxation parameter δ = 10 -1 . Various values of the time step are taken.

Figure 7 .

 7 Figure 7. Propagation of a solitary wave with relaxed schemes.

Figure 8 .

 8 Figure 8. Comparison of the exact and numerical solutions to the KdV using Sanz-Serna and extrapolated Relaxed and IMEX schemes.

Figure 9 .

 9 Figure 9. Comparison of the exact and numerical solutions to the KdV using Sanz-Serna and extrapolated Relaxed and IMEX schemes.

* In the case of continuous Galerkin methods this formulation is also sometimes called the modified one. † Here we mean I. M. Gelfand to avoid any confusion. J.-P. Chehab & D. Dutykh

Acknowledgments D. Dutykh would like to acknowledge the hospitality of the Laboratory LAMFA and of the University of Picardie Jules Verne during his visit in November 2016. Reciprocally, J.-P. Chehab acknowledges the hospitality of the Laboratory of Mathematics (LAMA UMR #5127) and of the University Savoie Mont Blanc during his visit in December 2017.

A. Other types of dispersive equations

In this Appendix we provide two other examples of relaxation for some widely used dispersive wave equations (scalar and system case).

A.1. Benjamin-Bona-Mahony equation

The celebrated Benjamin-Bona-Mahony (BBM) equation was derived first in [3,[START_REF] Peregrine | Calculations of the development of an undular bore[END_REF] and can be recast in the following dimensionless form:

Similarly to the KdV case the variable u (x, t) can be the free surface elevation or a horizontal fluid velocity either. It is possible to propose a similar relaxed formulation for the BBM equation as well. Indeed, let us recast it first in a conservative form:

The 2 nd derivative can be lowered by introducing additional variables:

to simplify the first equation:

The last step consists in adding relaxation terms to obtain an evolutionary system in all variables:

with δ ≪ 1 being again a small parameter.

A.2. Peregrine system

The Peregrine system was proposed by D. H. Peregrine (1967) in [START_REF] Peregrine | Long waves on a beach[END_REF]. In a particular case of a fluid layer of constant depth (i.e. the even bottom case) this system reads

where d > 0 is the constant fluid depth and g > 0 is the gravity acceleration. The variables η (x, t) and u (x, t) are the free surface elevation and the depth-averaged horizontal velocity correspondingly. The sketch of the fluid domain is shown in Figure 10. As the first step, the Peregrine system (A.1), (A.2) is rewritten in the following conservative form:

The last conservative Peregrine system can be lowered and relaxed similar to the BBM case (see the Section A.1) by introducing new variables:

where δ ≪ 1 and ρ (x, t) def := u (x, t) -1 3 d 2 w (x, t) is an auxiliary variable. If one takes the limit δ → 0 , we shall recover a formulation similar to what is used in local (or modified) continuous Galerkin methods [START_REF] Walkley | A finite element method for the one-dimensional extended Boussinesq equations[END_REF][START_REF] Walkley | A finite element method for the two-dimensional extended Boussinesq equations[END_REF].

Remark 2. The celebrated KdV and BBM models can be written for the variable u (x, t) being the free surface elevation η (x, t) or the horizontal velocity (after a few changes of variables). In this way these equations appear in the theory of water waves. However, the KdV and BBM-type equations appear in many other physical settings as well (see e.g. [11,17,[START_REF] Grimshaw | Internal Solitary Waves[END_REF]20,30]).