List of Figures

and for the Ozawa solution in Fig. 5.8 with N x = N y = 2 12 . . . . . . . 5.11 Evolution of max(|u| 2 ) and the numerically computed energy in dependence of time for a solution to the focusing DS II equation (3.2) for an initial condition of the form u(x, y, -6) = 1.1u l . . . . . . . . . . 5.12 Fourier coefficients at t = 0 for a solution to the focusing DS II equation (3.2) for an initial condition of the form u(x, y, -6) = 1.1u l . 5.13 Fourier coefficients at t = 6 for a solution to the focusing DS II equation (3.2) for an initial condition of the form u(x, y, -6) = 1.1u l . 5.14 Evolution of max(|u| 2 ) and the numerically computed energy in dependence of time for a solution to the focusing DS II equation (3.2) for an initial condition of the form u(x, y, -6) = 0.9u l . . . . . . . . . .

5.15

Solution to the focusing DS II equation (3.2) for an initial condition of the form u(x, y, -6) = 0.9u l for t = -3 and t = 0 in the first row and t = 3 and t = 6 below. . . . . . . . . . . . . . . . . . . . . . . . . 5. [START_REF] Besse | Numerical Study of the Davey-Stewartson System[END_REF] The Fourier coefficients at t = 0 of the solution to the focusing DS II equation (3.2) for an initial condition of the form u(x, y, -6) = 0.9u l . 5.17 an initial condition of the form (5.9) with κ = 0.9 for different times. 5.22 Evolution of max(|u| 2 ) and the numerically computed energy in dependence of time for the focusing DS II equation (3.2) for an initial condition of the form (5.9) with κ = 0.9. . . . . . . . . . . . . . . . . 5.23 Fourier coefficients of the solution to the focusing DS II equation (3.2) for an initial condition of the form (5.9) with κ = 0.9 at t = 0. . . . . 5. [START_REF] Certaine | The Solution of ordinary differential Equations with large Time Constants[END_REF] Contour plot for a solution to the focusing DS II equation (3.2) for an initial condition of the form (5.9) with κ = 1.1 for different times. for an initial condition of the form (5.12) with ν = 0.9. . . . . . . . . 5.38 Fourier coefficients of the solution to the focusing DS II equation (3.2) for an initial condition of the form (5.12) with ν = 0.9 at t = 0. . . . 5.39 Evolution of max(|u| 2 ) and the numerically computed energy E for a solution to the focusing DS II equation (3.2) for an initial condition of the form (5.12) with ν = 1.1. . . . . . . . . . . . . . . . . . . . . . 5.40 Fourier coefficients of the solution to the focusing DS II equation (3.2) for an initial condition of the form (5.12) with ν = 1.1 at t = 0.15. . . an initial condition of the form u(x, y, -3) = 1.1u l , DCRK method. . 6.4 Behavior of the Fourier Coefficients of u at t = 0 for an initial condition of the form u(x, y, -3) = 1.1u l , DCRK method. . . . . . . . . . . 6.5 Time evolution of max(|u| 2 ) and of the computed energy, for an initial condition of the form u(x, y, -3) = 0.9u l , DCRK method. . . . . . . . 6.6 Behavior of the Fourier Coefficients of u at t = 6 for an initial condition of the form u(x, y, -3) = 0.9u l , DCRK method. . . . . . . . . . . 6.7 Time evolution of the maximum of |u| 2 for an initial condition of the form (6.1), for A = 0.1, 0.

Abstract

Numerical analysis becomes a powerful resource in the study of partial differential equations (PDEs), allowing to illustrate existing theorems and find conjectures. By using sophisticated methods, questions which seem inaccessible before, like rapid oscillations or blow-up of solutions can be addressed in an approached way. Rapid oscillations in solutions are observed in dispersive PDEs without dissipation where solutions of the corresponding PDEs without dispersion present shocks. To solve numerically these oscillations, the use of efficient methods without using artificial numerical dissipation is necessary, in particular in the study of PDEs in some dimensions, done in this work. As studied PDEs in this context are typically stiff, efficient integration in time is the main problem. An analysis of exponential and symplectic integrators allowed to select and find the more efficient method for each PDE studied. The use of parallel computing permitted to address numerically questions of stability and blow-up in the Davey-Stewartson equation, in both stiff and non-stiff regimes.

L'analyse numérique se développe en un outil puissant dans l'étude des équations aux dérivées partielles (EDPs), permettant d'illustrer des théorèmes existants et de trouver des conjectures. En utilisant des techniques sophistiquées, des questions apparaissant inaccessibles avant, comme des oscillations rapides ou un blowup des solutions, peuvent être étudiées. Des oscillations rapides dans les solutions sont observées dans des EDPs dispersives sans dissipation où les solutions des EDPs correspondantes sans dispersion ont des chocs. Pour résoudre numériquement ces oscillations, l'application de méthodes efficaces introduisant peu de dissipation numérique artificielle est impérative, en particulier pour l'étude d' EDPs en plusieurs dimensions. Comme les EDPs étudiées dans ce contexte sont typiquement raides, l'intégration efficace dans le temps représente le principal problème. Une analyse des intégrants exponentiels et symplectiques a permi de déterminer les méthodes les plus efficaces pour chaque EDP étudiée. L'apprentissage et l'utilisation de techniques de parallélisation de codes numériques permet de nos jours de grandes avancées, plus précisément dans ce travail d'étudier numériquement la stabilité des solutions et l'apparition de blow-up dans l'équation de Davey-Stewartson. particulièrement les rapporteurs, le Pr. Joerg Frauendiener, et le Pr. Jean-Claude Saut. J'ai beaucoup apprécié également toutes leur suggestions de travail et je compte bien m'y attarder. D'un point de vue plus personnel, je remercie la famille Leslé, qui m'a acceuilli chez elle lors des evènements neigeux de décembre 2010 à Londres; ainsi que tous mes amis qui m'ont soutenu moralement, en particulier Sylvie, Nicolas et Thibault. Ma mère et mes beaux-parents m' ont aussi beaucoup encouragée, MERCI.

Enfin, à quoi rimerait tout ca si je ne remerciai pas celui sans qui, d'une part je n'aurai jamais fait une thèse, mais aussi sans qui je ne serai jamais arrivée au bout, j'aurai craqué avant... Mon mari, Jérôme Roidot, qui devrait recevoir une médaille pour m'avoir supportée durant ces trois années.

Introduction Physical and Mathematical Motivations

In every day life, wave phenomena are mainly known in the context of hydrodynamics. Consider for example a pebble falling into a lake in the below figure. One can see that concentric circles appear on the surface, they sprawl, and the place where the pebble hit the water becomes again rapidly undisturbed. The circular waves broaden whilst propagating in the medium (here a pond of finite depth), and the corresponding medium is said to be dispersive. In practice the waves are also damped due to dissipation. In this work we will, however, concentrate on the former effects and study nonlinear dispersive equations.

These arise in many fields of mathematical physics where wave phenomena play a role. Consequently this work is related to different domains of science as applied mathematics, computer sciences, fluid dynamics, general relativity, quantum mechanics... Equations of fluid dynamics have a prominent place in the theory of partial differential equations (PDEs). First because historically the dynamics of fluids could be easily observed, and secondly because the vast collection of complicated phenomena described by them as wave breaking and turbulence is modeled by nonlinear PDEs. The mathematical properties of such equations as singularities, instabilities,. . . appear also in other physical contexts described by similar nonlinear equations which gives PDEs from hydrodynamics an importance far beyond this field. Many mathematical properties of these equations are also present in more general nonlinear PDEs.

Among the most important of these PDEs are the Euler equations, which model in particular the classical problem of water waves.

Solutions to this nonlinear system of equations have complex properties which are yet under investigation. The nonlinearity in this system leads to the appearance of discontinuities called shock waves, as it can be seen by considering one of the simplest reduction of the Euler equations, the Hopf equation. Its solutions describe the distribution of velocities of a one-dimensional fluid without external forces, obviously a very simplified model. Even for regular initial data, solutions to the Hopf equation can develop singularities in finite time, i.e., the solutions may lead to a wave-breaking when the spatial derivatives blow up. The formal solution then becomes multi-valued, which is physically meaningless. This implies that dispersive effects neglected in the Hopf approximation which are given by terms with higher order spatial derivatives must be taken into account, which is achieved by considering the Korteweg-de Vries equation,

∂ t u + 6u∂ x u + 2 ∂ 3 x u = 0, x ∈ R, t ∈ R, (1)
which can also be derived from the Euler equations in certain regimes, and which models the propagation of one-dimensional waves in the limit of long wave lengths. These small dispersion effects lead to solutions which develop a zone of rapid modulated oscillations, called dispersive shocks and present in nature for instance as undular bores, in the region where the corresponding dispersionless equation has shocks. A typical example for the KdV equation is shown in the figure below, together with the physical observation of a related phenomenon in a river. Although their formation appears fascinating, these oscillations are mathematically difficult to describe, and represent an obstacle to the efficiency of numerical and even analytical methods.

Another interesting equation in this context is the nonlinear Schrödinger equation (NLS), given by

i ∂ t Ψ = -2 ∆ + λ|Ψ| p-1 Ψ, (2)
for Ψ a complex-valued function, with λ = ±1 and t ∈ R.

The latter is related to the well known Schrödinger equation in quantum mechanics,

i ∂ t Ψ = - 2 2 ∆ + V Ψ, (3)
in which > 0 is the Planck constant, ∆ is the Laplacian on R n , V (x) is a potential, and Ψ is a time dependent function in L 2 (R n , dx). It has a fundamental place in quantum mechanics 1 because it deals with the wave function giving a microscopic description of an object, for instance an electron. The transition to classical physics is formally achieved in the limit → 0, called the (semi)classical limit, and the parameter in (2) corresponds to the Planck constant in (3).

This limit is mathematically challenging since rapidly oscillating solutions of the Schrödinger equation appear which are also numerically problematic. But the semiclassical techniques provide exactly the mathematical tools to study the small dispersion limit of KdV and NLS equations (i.e., → 0), just this time for a nonlinear equation.

Both these purely dispersive equations have applications in the theory of water waves, but due to the universal importance of the equations, they appear for all sorts of wave phenomenona in acoustics, plasma physics, nonlinear optics... The semiclassical or small dispersion limit has been the subject of research in the last 40 years. Although an asymptotic description of these dispersive shocks is well known for certain integrable PDEs as KdV [START_REF] Lax | The small Dispersion Limit of the Korteweg de Vries Equation I, II, III[END_REF][START_REF] Venakides | The Korteweg-de Vries Equation with small Dispersion: higher order Lax-Levermore Theory[END_REF][START_REF] Deift | New result in small dispersion KdV by an extension of the steepest descent method for Riemann-Hilbert problems[END_REF] and the NLS equation for certain classes of initial data [START_REF] Jin | The behavior of solutions of the NLS equation in the semiclassical limit[END_REF][START_REF] Kamvissis | Semiclassical Soliton Ensembles for the Focusing Nonlinear Schrödinger Equation[END_REF][START_REF] Tovbis | On semiclassical (zero Dispersion Limit) Solutions of the focusing nonlinear Schrödinger Equation[END_REF], no such description is known for (2 + 1)-dimensional PDEs which, in addition, can have solutions which blow up, i.e., solutions have after finite time lower regularity than the initial data. It is known for many of these PDE when blowup can occur, but for the precise mechanism of the blow-up not even conjectures exist.

Even though most of the research so far is theoretical, there is a lack of results in the study of such phenomena, in particular in several space dimensions. At this point, exchanges between theory and numerical simulations are necessary. It is important to point out that a constructive theory can more easily be formulated with the use of mathematical modeling on computers and computational experiments. This allows to check existing theories, and also to give conjectures which may serve to define new directions in the development of new theories. The advances in the construction of powerful parallel computers (clusters), allow nowadays to investigate large and complex problems which could not be addressed before.

Special attention in this work is given to (2+1)-generalizations of the NLS and the KdV equations, which are completely integrable, the Kadomtsev-Petviashvili equation, and the Davey-Stewartson equations.

The Kadomtsev-Petviashvili equation reads

∂ x ∂ t u + 6u∂ x u + 2 ∂ xxx u + λ∂ yy u = 0, λ = ±1 (4)
where (x, y, t) ∈ R x × R y × R t and where 1 is a small scaling parameter. The case λ = -1 corresponds to the KP I model and the case λ = 1 to the KP II model. KP I is used to model water waves when the surface tension is strong, and KP II when the surface tension is weak. KP equations also appear in other physical fields, for example in ferromagnetic media [START_REF] Turitsyn | Stability of magneto-elastic solitons and selffocusing of sound in antiferromagnets[END_REF] to model sound waves and in Bose-Einstein condensates [START_REF] Huang | Two-dimensional solitons in Bose-Einstein condensates with a disk-shaped trap[END_REF] to describe nonlinear matter-wave pulses.

The Davey-Stewartson equations are given by

i∂ t u + ∂ xx u -α∂ yy u + 2ρ φ + |u| 2 u = 0 ( 5
)
∂ xx φ + α∂ yy φ + 2∂ xx |u| 2 = 0 (6)
for the (complex) amplitude u(x, y, t) and for the (real) mean velocity potential φ, with α, ρ = ±1. Equation ( 5) is a nonlinear Schrödinger equation with a forcing term, and equation ( 6) is a linear equation which is either elliptic or hyperbolic, depending on the sign of α. Thus, these systems have been classified in [START_REF] Ghidaglia | On the initial value problem for the Davey-Stewartson systems[END_REF], with respect to the second order differential operators in equations ( 5) and ( 6). The hyperbolic-elliptic case (α = 1) is known as DS II equation, and the DS I equation corresponds to the elliptic-hyperbolic case (α = -1). Both cases arise as a higher dimensional generalization of the cubic NLS equation (p = 3 in ( 2)), and also from physical considerations. For example they model the evolution of weakly nonlinear water waves that travel predominantly in one direction, but in which the wave amplitude is modulated slowly in two horizontal directions [START_REF] Davey | On three-dimensional packets of surface waves[END_REF], [START_REF] Djordjevic | On two-dimensional Packets of Capillarity-Gravity Waves[END_REF]. They are also used in plasma physics [START_REF] Nishinari | A new-type of Soliton Behavior of the Davey-Stewartson Equations in a Plasma System[END_REF][START_REF] Nishinari | Multi-dimensional localized Behavior of electrostatic Ion Wave in a magnetized Plasma[END_REF], to describe the evolution of a plasma under the action of a magnetic field. The focusing cubic NLS equation (λ = -1) has the critical nonlinearity in 2 space dimensions with respect to finite time blow-up. A question of particular interest is thus to explore the blow-up behavior of DS solutions, and its relation to the NLS case.

In view of the importance of the equations, and the difficulty of the open mathematical questions, efficient numerical algorithms are needed to enable extensive numerical studies of the PDE. Since critical phenomena are generally believed to be independent of the chosen boundary conditions, we study a periodic setting. Such settings also include rapidly decreasing functions which can be periodically continued within the finite numerical precision. This allows to approximate the spatial dependence via truncated Fourier series which leads for the studied equations to large systems of ODEs. The use of Fourier methods not only gives spectral accuracy in the spatial coordinates, but also minimizes the introduction of numerical dissipation which is important in the study of dispersive effects. In addition, they usually require less time and memory than other alternative approaches, such as finite differences or finite elements methods. This 'memory-minimizing' property is crucial in the implementation of parallel codes, because the memory used by an algorithm has a first-order effect on its runtime (time of completion) on modern parallel systems. The fact that serial codes to compute discrete Fourier transforms are already implemented and obtain high-performance (since the well known fast Fourier transformation (FFT) algorithm of Cooley and Tukey (1965)) also makes the choice of spectral methods convenient, since the FFT can be directly used for the code.

In Fourier space, equations ( 4) and ( 5)-( 6) have the form

v t = Lv + N(v, t), ( 7
)
where v denotes the (discrete) Fourier transform of u, and where L and N denote linear and nonlinear operators, respectively. The resulting systems of ODEs are classical examples of stiff equations where the stiffness is related to the linear part L due to the high order of derivatives, whereas the nonlinear part contains only low order derivatives. In the small dispersion limit, this stiffness is still present despite the small term 2 in L. This is due to the fact that the smaller , the higher wavenumbers are needed to resolve the rapid oscillations. A system can be called stiff if it involves a wide range of timescales. For example, the solution of the following system,

ẋ = M x, with M = -1 0 0 -µ (8)
with µ positive and large, is x(t) = (x 1 (0)e -t , x 2 (0)e -µt ), for which the second component tends to zero faster than the first one. These separated time scales in the system characterize stiff systems. The latter are difficult to solve, because the suppression of exponentially growing modes (stability of the used scheme) implies that the time steps have to be chosen prohibitively small for conventional explicit numerical methods.

Several approaches to deal efficiently with equations of the form (2.17) with a linear stiff part can be found in the literature: implicit-explicit (IMEX) methods, in which the main idea is to solve the linear part (responsible for the stiffness) by an implicit method, and the nonlinear one by an explicit method; time splitting methods, in which the initial equation is split into two or more equations explicitly integrable; integrating factor (IF) and deferred correction schemes as well as sliders and exponential time differencing methods may also represent an efficient tool to this kind of equations. These methods are explained in detail in chapter 2. To avoid as much as possible a pollution of the Fourier coefficients by errors due to the finite difference schemes for the time integration and to allow for larger time steps, we only consider fourth order schemes as promising candidate methods for our purposes. While standard explicit schemes would impose prohibitively small time steps due to stability requirements, stable implicit schemes would be computationally too expensive in 2 + 1 dimensions. We therefore only consider explicit schemes, and select the ones we will use in the second chapter.

The KdV equation was the first nonlinear evolution equation discovered to be integrable by the inverse scattering approach (the second being the NLS equation), and their generalizations under consideration in this work are also known to be 'completely integrable, or 'exactly solvable', in a sense to be discussed in the first chapter. This complete integrability allows Hamiltonian formulations and the construction of exact solutions, typically solitons, which thus provide popular test cases for numerical algorithms. But as we will show the latter often test the equations in a non-stiff regime. The main challenge in the study of critical phenomena as dispersive shocks and blowup is, however, the numerical resolution of strong gradients in the presence of which the above equations are stiff. This implies that algorithms that perform well for solitons might not be efficient in the context studied here. An identification of the most efficient method for each equation is thus imperative, and discussed in the third chapter.

For the focusing DS II equation, it is crucial to provide sufficient spatial resolution for the peaks. The reason is the modulational instability of this equation, that is self-induced amplitude modulation of a continuous wave propagating in a nonlinear medium. This instability leads to numerical problems if there is no sufficient spatial resolution of the maxima, since the high Fourier modes will grow in this case and eventually destroy the numerical solution. In [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger Equation[END_REF], a resolution of 2 13 modes was necessary for the study of the semiclassical limit of the focusing NLS in 1 + 1 dimensions. The possibility of blowup in DS requires at least the same resolution.

Despite the perpetual advances in the construction of more and more powerful computers, an essential tool nowadays for numerical simulations is the parallel computing. It allows to address more complex and larger problems than before, and constitutes, in our opinion, the tool that research needs to advance, in particular in cases where theories do not yet succeed. We thus wrote parallelized codes in Fortran 90 to perform a valid numerical study of this equation.

Outline of the thesis

The first chapter is devoted to the presentation of dispersive PDEs and related physical phenomena of interest in this work, as dispersive shock waves and blow-up. Some aspects of the complete integrability of a PDE are reviewed, as well as some analytical aspects and explicit solutions of the KP and DS II equation.

Chapter II gives an introduction to stiff systems and numerical methods used in the later numerical experiments. We first select a method to discretize the differential equation in space and incorporate boundary conditions, before integrating equations in time by using an appropriate method chosen in accordance with the properties of the spectrum of the discrete operator resulting from the discretization in space. We thus discuss some candidate methods to deal with stiffness and dispersive shocks.

In chapter III, we compare the performance of several time integration schemes to deal with the KP and DS II equations, and identify in each case the most efficient method. It will be shown that fourth order time stepping schemes can be efficiently used here. As expected, algorithms that perform well for solitons are in general not efficient in the small dispersion limit. We found that in the non-stiff regime, Driscoll's method [START_REF] Driscoll | A composite Runge-Kutta Method for the spectral Solution of semilinear PDEs[END_REF] performs well, whereas in the stiff regime, Exponential Time Differencing methods are more efficient. We observe numerically the phenomenon of order reduction found by Hochbruck and Ostermann in [START_REF] Hochbruck | Exponential Runge-Kutta Methods for semilinear parabolic Problems[END_REF] for parabolic systems, and establish that the numerical conservation of the L 2 norm is a good indicator for the quality of numerics. Implicit schemes as a two-stage Gauss scheme are computationally too expensive in the stiff regime, whereas standard explicit schemes as Runge-Kutta require for stability conditions prohibitively small time steps for the KP and DS equations. In the small dispersion limit, the KP equation shows a similar behavior as the KdV equation in the one-dimensional case. We observe dispersive shocks where the corresponding dispersionless equation has shocks.

The modulational instability of the DS equation requires high spatial resolution we could not achieve in Matlab on the used computers for general initial data in chapter III. Thus, in the (short) chapter IV, two codes have been parallelized to allow the use of higher spatial resolution without allocating too much memory per processor. We give there basic notions of parallel computing, and explain how we develop our parallel codes.

In chapter V, we use a parallel code for the time-splitting scheme of chapter III for the focusing DS II equation. We explore the formation of blow up in DS II, and study perturbations of exact solutions of this equation. We illustrate that splitting codes produce some artificial continuation of the solution beyond the blow-up time. The latter can be identified by a jump of the energy, a conserved quantity of DS. We found, in accordance with the conjecture in [START_REF] Mcconnell | Localised coherent solutions of the DSI and DSII equations-a numerical study[END_REF], that solutions to the focusing DS II equation either blow up or disperse. In particular the lump is unstable against both blowup and dispersion, in contrast to the lump of the KP I equation that appears to be stable.

In chapter VI, the comparison of the time splitting scheme with the most efficient method selected in the chapter III for the study of the DS II equation, the method presented in [START_REF] Driscoll | A composite Runge-Kutta Method for the spectral Solution of semilinear PDEs[END_REF] confirms the conclusion of chapter III by showing higher efficiency than the splitting scheme. Finally we investigate the semiclassical limit of this system. We find that the focusing DS II equation in the semiclassical limit leads to a blow up as in the critical one-dimensional case (in which the nonlinearity is quintic), for initial data with radial symmetry, whereas more general initial data lead to solutions with dispersive shocks.

Finally, we do a conclusion and present future axis of research.

Chapter 1

Dispersive Partial Differential Equations and Integrable Systems

We introduce in this chapter the notion of dispersive partial differential equations (PDEs), both in the linear and the nonlinear case. We then discuss complete integrability of PDEs and related properties. After a presentation of the physical effects due to the nonlinearity of dispersive PDEs under consideration in this work, we focus on two (2+1)-dimensional equations for which we review some analytical facts and explicit solutions which will be of importance in the following.

Linear Dispersive Partial Differential Equations

We introduce in this section the notion of dispersive PDEs for linear equations.

In the following, t represents the time variable, and the components of x in R d are spatial variables. We consider linear evolution equations with constant coefficients, for a function (or even a distribution)

u : R + × R d → V (t, x) → u(t, x) (1.1)
where V is a finite-dimensional Hilbert space. We mainly illustrate in this part first order PDEs in t, having as spatial domain R d or T d = (R/2πZ) d , i.e., PDEs of the form

∂ t u = Lu, u(0, x) = u 0 (x) (1.2)
where L is a skew adjoint differential operator in space with constant coefficients, of the form

Lu(x) := |α|≤m c α ∂ α x u(x) = 0 (1.3)
where m ∈ N (which represents the order of the differential operator),

α := (α 1 , ..., α d ) ∈ Z d
+ ranges over all multi-indices with |α| := d j=1 α j (assumed to be bounded by m), c α are coefficients (real or complex) that do not depend on x, and ∂ α x is the partial derivative of order |α| defined by

∂ α x := ∂ α 1 x 1 ...∂ α d x d . (1.4)
A plane wave is a wave which only depends on t and ξ• x, for t ξ ∈ R d , called the wave vector; the symbol ( • ) represents here the scalar product on L 2 (R d ) defined by

u• v = R d uv dx (1.5)
By definition, a plane wave is constant on any planes orthogonal to the wave vector.

A traveling plane wave is a plane wave of the form

u : (t, x) → u(t, x) = U (ξ• x -ωt), (1.6)
where ω ∈ R and U is constant. Its velocity is c := ω |ξ| , which means that

u(t, x) = u(0, x -ctn), with n := ξ |ξ| .
A monochromatic plane wave is a particular plane wave periodic in the phase of the form u :

(t, x) → u(t, x) = U e i(ξ•x-ωt) (1.7)
The quantity ω is called the phase of the wave, ω 2π is its frequency (i.e., the inverse of the temporal period), and λ = 2π |ξ| its wave length (i.e., its spatial period in the direction n).

By looking for monochromatic plane wave (MPW) solutions of a linear PDE, one gets an equation for (ω, ξ) called the dispersion relation.

When the latter is solvable in the form ω = ω(ξ), one can define two velocities:

The phase velocity,

v p := ω(ξ) |ξ| ξ |ξ| ,
and the group velocity,

v g := ∇ω,
where ∇ is acting on ξ.

When v p = v g , the PDE is called dispersive, and its solutions are dispersive waves. The PDE is also called dispersive when v p depends on ξ.

To illustrate this notion, we consider the following well known (linear) equations:

1-d Transport equation:

∂ t u + αu = 0, α ∈ R Wave equation: ∂ 2 tt u -α 2 ∆u = 0, α ∈ R Klein-Gordon equation: ∂ 2 tt u -α 2 ∆u + βu = 0, α ∈ R + , β ∈ R * Airy equation: ∂ t u + α∂ x u + β∂ 3 x u = 0, α, β ∈ R, Schrödinger equation: i∂ t u + ∆u = 0, > 0
The 1-dimensional transport equation has (complex) MPW solutions of the form u(t, x) = U m e i(ξx-ωt) if ω is linked to ξ by ω -αξ = 0 (dispersion relation), leading to ω = αξ.

Their phase velocity being constant, v p = α, the 1-dimensional transport equation is not dispersive.

The dispersion relation of the wave equation is ω 2 -α 2 |ξ| 2 = 0, which implies ω = ±α|ξ|, and

v p = ±α ξ |ξ| = v g
Thus, in the one-dimensional case, the wave equation is not dispersive, we have v p = v g ∈ {-α, α}, and wave train (both on the right and on the left) are not dispersive, see the illustration in the Fig. 1.1 for α = 1. However, for higher dimensional cases, the phase velocity then depending on ξ, the wave equation is weakly dispersive.

For the Klein-Gordon equation, which is the wave equation plus a damping term (βu, β > 0), one finds

ω = ± β + α 2 |ξ| 2 , v p = ± β + α 2 |ξ| 2 |ξ| ξ |ξ| and v g = ± α 2 ξ β + α 2 |ξ| 2 .
The Klein-Gordon equation is thus dispersive.

The Airy equation (also referred to the linear Korteweg-de Vries equation) is one of the main examples of linear dispersive equations, for which we have ω = αξ -βξ 3 , v p = α -βξ 2 , and v g = α -3βξ 2 ; together with the Schrödinger equation, for which one gets

ω = |ξ| 2 , v p = ± |ξ|, and v g = 2 |ξ| , see Fig. 1.2.

Remark 1.1

The above description is specific to linear PDEs with constant coefficients. Another approach to present these notions, and allowing the study of linear (and even nonlinear) PDEs with variable coefficients, is based on geometric optic. The idea is to look for solutions to the PDE depending on a small parameter , written in dependence of the slow variables x → x and t → t. In fact, one looks for solutions in the form of wave trains, with an amplitude of the order of O(1) and a phase of order of O(1/ ), i.e., u(x, t) = A( x, t; )e iS( x, t)/ .

Linear PDEs are generally approximations of nonlinear models that are much more difficult to solve, and we present in the next section a few semilinear dispersive PDEs of great importance in mathematical physics.

Semilinear Dispersive PDEs

In this section, we focus on nonlinear dispersive equations, which mainly describe wave propagation in a weakly nonlinear and dispersive environment.

Let Ω be a domain of R n with boundary Γ. We denote by u = u(t, x) a function from R × Ω to R p . Typically, we study equations of the form

∂ t u = Lu + N [u] (1.8)
where N represents the nonlinear part, which may contain low order spatial derivatives of u, and where the linear operator L is defined as before. Such PDEs in which the deriva-tives in the nonlinearity are of lower order than in the linear term are called semi-linear.

We present here several typical examples of semilinear dispersive PDEs which appear in many domains of physics and have been intensively studied for many years. We split them into two main families: equations of Nonlinear Schrödinger (NLS) type and equations of Korteweg-de Vries (KdV) type.

Equations of NLS Type.

The Cauchy problem associated to the Nonlinear Schrödinger equation (NLS) is given by

i∂ t u + ∆u = λ |u| p-1 u, u(0, x) = u 0 (x) (1.9)
where u(t, x) is a complex-valued function, λ = ±1, t ∈ R, x ∈ R d and the initial data u 0 is specified and lies in a given Sobolev space

H s x (R d ) (or H s-1 x (R d ))
. The exponent 1 < p < ∞ represents the power of the nonlinearity. When p is an odd integer (algebraic cases), the nonlinearity z ∈ C → |z| p-1 z is smooth, and naturally associated to a Hamiltonian potential V (u) := λ p+1 |u| p+1 . The cases p = 3 and p = 5 are referred to as the cubic case and the quintic case respectively, which are the most common and studied cases in physical models.

The parameter λ denotes whether the nonlinearity is defocusing, or focusing (sometimes referred to as repulsive and attractive respectively). More precisely, this terminology can be first justified by looking at explicit solutions obtained by the method of separation of variables.

The plane wave u(t, x) = e iξ•x v(t), with |v(t)| = α is a solution of the NLS equation if v(t) satisfies the ODE

∂ t v = -i(|ξ| 2 + λ|v| p-1 )v,
which is explicitly solvable, leading to

u(t, x) = αe iξ•x e i|ξ| 2 t e iλ|α| p-1 t
for any α ∈ C and ξ ∈ Z d . Both terms e i|ξ| 2 t and e iλ|α| p-1 t contribute to the oscillations in time. Thus, in the defocusing case (λ = 1), both oscillations are anti-clockwise, leading to an amplification of the dispersive effect of the linear equation. In the focusing case (λ = -1) the corresponding term is instead trying to suppress the latter. If the amplitude α is small compared to the frequency ξ, the dispersive effect is stronger, and when α is large, the focusing effect takes over.

The NLS equation is a purely dispersive equation with important applications in nonlinear optics, hydrodynamics, and plasma physics (where (1.9) is viewed as a particular case of the Zakharov equation).

Another example of a dispersive nonlinear PDE of NLS type is the integrable Davey Stewartson system (DS), given by:

i∂ t u + ∂ xx u -α∂ yy u + 2ρ φ + |u| 2 u = 0 (1.10) ∂ xx φ + α∂ yy φ + 2∂ xx |u| 2 = 0 (1.11)
for the (complex) amplitude u(x, y, t) and for the (real) mean velocity potential φ, with α, ρ = ±1. Equation (1.10) is a nonlinear Schrödinger equation with a forcing term, and equation (1.11) is a linear equation which is either elliptic or hyperbolic, depending on the sign of α. Thus, these systems have been classified in [START_REF] Ghidaglia | On the initial value problem for the Davey-Stewartson systems[END_REF], with respect to the second order differential operators in equations (1.10) and (1.11). The hyperbolic-elliptic case (α = 1) is known as DS II equation, and the DS I equation corresponds to the elliptic-hyperbolic case (α = -1). Both cases arise as a higher dimensional generalization of the nonlinear Schrödinger equation (NLS), and also from physical considerations, for example they model the evolution of weakly nonlinear water waves that travel predominantly in one direction, but in which the wave amplitude is modulated slowly in two horizontal directions [START_REF] Davey | On three-dimensional packets of surface waves[END_REF], [START_REF] Djordjevic | On two-dimensional Packets of Capillarity-Gravity Waves[END_REF]. They are also used in plasma physics [START_REF] Nishinari | A new-type of Soliton Behavior of the Davey-Stewartson Equations in a Plasma System[END_REF][START_REF] Nishinari | Multi-dimensional localized Behavior of electrostatic Ion Wave in a magnetized Plasma[END_REF], to describe the evolution of a plasma under the action of a magnetic field.

We study in more detail these systems in the following.

Equations of KdV Type.

These equations can be written in the form

∂ t u + ∂ x (M u) + ∂ x (f (u)) = 0 (1.12)
where x ∈ R n , t > 0, u(x, t) is a real valued function and M is defined as

F(M u)(ξ) = ς(ξ)F(u)(ξ),
with ς real-valued and F the n-dimensional Fourier transform.

The typical example is the Korteweg de Vries equation (KdV) given by

∂ t u + 6u∂ x u + ∂ 3 x u = 0, x ∈ R, t ∈ R (1.13)
which is a purely dispersive PDE describing one-dimensional wave phenomena in the limit of long wave lengths.

A second interesting dispersive nonlinear PDE of KdV type is the Benjamin-Ono equation (BO) given by

∂ t u + u∂ x u + H[∂ xx u] = 0, x ∈ R, t ∈ R (1.14)
where H is the Hilbert Transform Operator defined by 1 :

H[f (x)] = 1 π P.V. ∞ -∞ f (y) y -x dy (1.15)
The BO equation was introduced by Benjamin (1967) and Ono (1975) and describes internal waves of deep-stratified fluids.

A last example of dispersive nonlinear equation of KdV type is a generalization to 2space dimensions due to Kadomtsev and Petviashvili in [START_REF] Kadomtsev | On the stability of solitary waves in weakly dispersing media[END_REF], who studied the stability of the KdV soliton against weak transverse perturbations: the Kadomtsev-Petviashvili equation (KP) which reads in dimensionless form

∂ x ∂ t u + 6u∂ x u + ∂ 3 x u + 3λ∂ 2 y u = 0, (x, y) ∈ R 2 , t ∈ R (1.16)
with λ = ±1. The case λ = -1 corresponds to the KP I model and is known to have a focusing effect, whereas the case λ = 1 called the KP II model is defocusing. KP I is used to model water waves when the surface tension is strong, and KP II when the surface tension is weak. KP equations also appear in other physical fields, for example in ferromagnetic media [START_REF] Turitsyn | Stability of magneto-elastic solitons and selffocusing of sound in antiferromagnets[END_REF] to model sound waves and in Bose-Einstein condensates [START_REF] Huang | Two-dimensional solitons in Bose-Einstein condensates with a disk-shaped trap[END_REF] to describe nonlinear matter-wave pulses.

We will focus in this work essentially on the DS system and the KP equations, for which we will give some review on analytical properties in the following. All dispersive PDEs cited above are known to be completely integrable by Inverse Scattering Transform methods, which leads to interesting properties and to exact solutions for which the energy density is mainly localized in a small part of space, and which evolve without change of form (so-called soliton solutions), that we review in the next section.

Complete Integrability

The concept of completely integrable system arose in the 19th century in the context of finite-dimensional classical mechanics with the goal to find exact solutions to Newton's equations of motion.

Hamilton had reformulated Newton's equations by introducing the so-called canonical coordinates (q 1 , ..., q n ) (generalized positions) and (p 1 , ..., p n ) (generalized momenta) to describe a mechanical system with n degrees of freedom.

The temporal evolution of any initial state (x 0 , p 0 ) is then governed by Hamilton's equations of motion, qj = ∂H ∂p j ṗj = -∂H ∂q j j = 1...n, 1 P.V. stands for principal value defined in a region Ω of the space R 2n (q 1 , ..., q n , p 1 , ..., p n ), called the phase space, with the function H(q, p) ∈ C ∞ (Ω) representing the total mechanical energy of the system and called the Hamiltonian of the system.

We define the canonical Poisson bracket for smooth functions f and g on the phase space Ω by

{f, g} = n j=1 ∂f ∂q j ∂g ∂p j - ∂f ∂p j ∂g ∂q j ,
for which we have the canonical relations

{q i , q j } = 0, {p i , p j } = 0, {q i , p j } = δ i,j , i, j = 1, ..., n
where δ i,j denotes the Kronecker symbol.

Hamilton's equations can be rewritten as qj = {q j , H} ṗj = {p j , H} j = 1...n.

The most important objects in solving the equation are the first integrals of motion, also called conserved quantities.

A function f ∈ C ∞ (Ω) is a first integral of motion if {f, H} = 0.
A dynamical system defined by a given Hamiltonian H on a 2n-dimensional phase space Ω is called completely integrable if it has n independent integrals of motion in involution (f 1 , ..., f n ) on Ω (sometimes also referred to as Hamiltonians), i.e., all Poisson brackets {f j , f k }, j, k = 1..n vanish. Thus these first integrals are conserved under the Hamiltonian evolution on Ω generated by each of them.

There is a key result, the Liouville-Arnold theorem, which ensures that for completely integrable systems, there exists a canonical transformation to action-angle coordinates, a special set of variables (J i , φ i ) on the phase space such that the transformed Hamiltonians depend only on the action variables J, which are a combination of integrals of motion, and the motion is constrained to the surface of a torus, known as the invariant Liouville torus.

The coordinates on the torus are the angle variables φ i , see [START_REF] Arnol | Dynamical Systems[END_REF]. Consequently, the actions are conserved, and the angles evolve linearly in the evolution parameters (t 1 , ..., t n ), contrary to (q, p) for which the dependence in t is highly nonlinear for typical systems. Hence one can explicitly solve Hamilton's equations in such cases.

From the Finite to the Infinite Dimensional Notion of Complete Integrability

After a lack of interest in the subject, most probably due to the results by Poincaré, which stated that for the systems usually considered in classical mechanics (Hamiltonian systems), integrability is a highly exceptional property, the discovery of the soliton phenomena by Zabusky and Kruskal (1965, see [START_REF] Zabusky | Interaction of Solitons in a collisionless Plasma and the Recurrence of initial States[END_REF]) led to a revival of the domain and to the emergence of new results and general structures in the domain.

Except for some insights concerning the KdV equation (1.13), soliton theory did not see much progress, although its empirical discovery dates back to 1834, when J.S. Russell first observed the solitary wave, a hump of water moving with constant speed and shape along a canal. The corresponding 1-soliton solution to the KdV equation

u(x, t) = c 2 sech 2 √ cξ 2 ,
where c > 0 is the wave speed and ξ = x -ct is a moving coordinate, also dates back to the 19th century.

Computer simulations by Zabusky and Kruskal, who studied collisions of n solitary waves led to the discovery of extraordinary stability properties of these solitary waves (see [START_REF] Zabusky | Interaction of Solitons in a collisionless Plasma and the Recurrence of initial States[END_REF]). They observed that these waves emerge with the same velocities and shapes as before the collision. This was surprising, since the KdV equation being nonlinear, solutions cannot be linearly superposed. Actually, the presence of a nonlinear interaction was understood with the explicit form of the solutions: the positions of the solitary waves are shifted, compared to the positions arising from a linear superposition. These particle-like properties led to the name of soliton for such waves.

The connection with the concept of completely integrable system was first made by Zakharov and Faddeev (see [START_REF] Zabusky | Korteweg-de Vries Equation, a completely integrable Hamiltonian System[END_REF]), who had shown that the KdV equation has an infinite number of conservation laws, and that there exists a linearizing transformation, which maps the initial value u(0, x) (which is assumed to decay sufficiently rapidly at infinity) for the KdV Cauchy problem, to spectral and scattering data of the Schrödinger operator

d 2 dx 2 -u(0, x).
The problem is thus transformed into a linear time evolution of these data, and the solution u(t, x) can be constructed via the inverse map, the so-called Inverse Scattering Transform (IST). Zakharov and Faddeev showed that the KdV equation may be viewed as an infinite-dimensional classical integrable system, the spectral and scattering data being the action-angle variables, the IST the (inverse of the) action-angle map, and the infinity of conserved quantities the Poisson commuting Hamiltonians.

Since then, the number of nonlinear partial differential equations which have been shown to have Hamiltonian form on appropriate infinite dimensional manifolds and to have an infinite number of conservation laws has increased.

The infinite dimensional generalization of the Hamiltonian formalism is based on two important functionals: the variational gradient, and the Poisson bracket. It is important to note that in contrast to the finite dimensional case, no concept of a canonical structure exists.

A PDE of the (abstract) form

∂u ∂t = P (x, u, du, ..., d m u), x ∈ R d (1. 17
)
where P is regular and d denotes the differentiation only with respect to x, admits an Hamiltonian structure if it can be written in the form

∂u ∂t = J δH[u], (1.18)
where

H[u] := R d H(x, u, du, ..., d m u)dx
is the Hamiltonian. Its variational gradient δH is characterized by

dH[u + θh] dθ | θ=0 = δH[u]• h dx, ( 1.19)
and can be expressed via the Euler operator by

δH :=    E 1 H . . . E n H    , (1.20)
where

E j H := α (-1) |α| ∂ α ∂H ∂ u j,α . (1.21) For example, if d = m = 1 and H = H(u, ∂ x u), then EH = ∂H ∂u -∂ x ∂H ∂ux .
Finally, J is the differential operator (eventually nonlinear) defined by

(J u) (x) := J(x, u, du, ..., d m u)
which has to satisfy the two following properties:

Anti-symmetry : For all pairs (Q, R) of functionals,

δQ[u]• J δR[u] = -δR[u]• J δQ[u]
for all regular functions u such that both sides of the equality are well defined. In other words, the bracket {, } defining the functional {Q, R} by

{Q, R}[u] := δQ[u]• J δR[u]
is anti-symmetric.

Jacobi Identity : For all triplets (Q, R, S) the bracket {, } satisfies {Q, {R, S}} + {R, {S, Q}} + {S, {Q, R}} = 0 Thus, {, } is a Poisson bracket.

Such equations are referred to completely integrable, or exactly solvable, although it no longer makes sense to count half the number of dimensions2 .

Soon afterwards, additional structural features have been shown to be present for many of these equations, including the existence of infinitely many conserved quantities, a Lax pair formulation, for which the connection to integrability was made by Faddeev, Zakharov and Gardner (see [START_REF] Zabusky | Korteweg-de Vries Equation, a completely integrable Hamiltonian System[END_REF]); Bäcklund transformations, a prolongation algebra, and a linearizing map playing a role comparable to Fourier transformation for linear PDEs.

The Lax pair representation of an equation , see [START_REF] Lax | Integrals of nonlinear equations of evolution and solitary waves[END_REF], [START_REF] Lax | Scattering theory[END_REF], consists in presenting the equations of motion of the system in the form

d dt A λ = [A λ , B λ ],
where entries of A λ and B λ are functions of the dynamical variables and depend also on λ, the spectral parameter, and where [, ] denotes the commutator of matrices. For example, the Lax representation of the KdV equation (1.13) is given by

∂ xx Ψ = (u -λ)Ψ (1.22) ∂ t Ψ = ∂ 3 x + 3 2 u∂ x Ψ + 3 4 ∂ x uΨ + αΨ (1.23)
where the solution u(t, x) of the KdV equation plays the role of a scattering potential, Ψ(t, x, λ) is the corresponding eigenfunction, λ is the spectral parameter and α is an arbitrary constant.

The Lax representation leads to first integrals of motion, and to the concept of a spectral curve.

There is no fully satisfactory definition of complete integrability for classical systems with infinitely many degrees of freedom, making it difficult to determine the integrability of a given system of equations. Instead, the term is used when some of the above structural features are present. Thus higher-dimensional PDEs with soliton solutions have been shown to be completely integrable, and in particular the Kadomtsev-Petviashvili and the Davey-Stewartson equations.

The special character of integrable PDEs also shows up in associated ODEs: the only movable singularities of the latter are poles (the Painlevé property), see [START_REF] Zakharov | What is Integrability?[END_REF].

Completely integrable systems are related to many fields in mathematics. For example the construction and study of the above-mentioned linearizing map leads to questions in functional analysis (scattering and inverse scattering theory, spectral theory, integral equations,...), function theory and algebraic geometry (Riemann-Hilbert problem, elliptic and hyperelliptic Riemann surfaces, Baker-Akhiezer and theta functions, Jacobi varieties,...), differential geometry (geodesic flows on groups and symmetric spaces), symplectic geometry (Hamiltonian structures, moment map, Weinstein-Marsden reductions,...), and in Lie algebra and Lie group theory.

From the perspective of dynamical systems, the existence of stable solitary wave solutions for certain completely integrable PDEs may be viewed as a consequence of a delicate balance between two competing effects: the nonlinearity of the equation, which tends to focus the waves, and the linear dispersive part which tends to smear them out.

We illustrate this now at the example of the generalized KdV equation.

Solitonic Solutions, an Equilibrium between Dispersion and Nonlinearity.

The generalized KdV equation (gKdV) is given by

∂ t u + ∂ xxx u + u p u x = 0,
where p is a positive integer.

The dispersive effect of the linear part can be seen from the linearized gKdV equation, the Airy equation,

∂ t u + ∂ xxx u = 0, u(0, x) = u 0 , (1.24)
which represents one of the simplest dispersive equation. Its solutions have the general form of a superposition of oscillatory waves with a dispersion relation ω = -ξ 3 . These waves have a phase velocity depending on ξ, which characterizes a dispersive medium, and thus u(t, x) describes dispersive waves whose profile develops an oscillatory decaying tail as we can see in Fig. 1.3, in which the initial condition for the Cauchy problem (1.24) is chosen as

u 0 = -sech 2 x.
The steepening effect of the nonlinear part can be seen by considering the initial value problem associated to the nonlinear part of the gKdV equation, which is a quasi-linear hyperbolic equation.

∂ t + u p ∂ x u = 0, u(0, x) = u 0 (x), p ∈ N + * , ( 1.25)
Solving by the characteristic method

dt dτ = 1, dx dτ = u p , du dτ = 0, with initial conditions t(0) = 0, x(0) = s and u(0, s) = f (s),
one gets for the general solution u(t, x) = f (x -u p t). This solution holds until it breaks down, after the collision of characteristics at

t c = min x∈R - 1 ∂ x f (x)
where both ∂ t u and ∂ x u become singular. There is a blow-up of the gradient, a gradient catastrophe in t c . There cannot exist a classical solution to (1.25) after t c .

As an example, we consider the equation (1.25) with the initial data t = 0, u 0 =sech 2 x and n = 1. Figure 1.4 shows the evolution of the solution in time. As time increases, the slope of the wave steepens to the right and flattens to the left. After the breakup time, the solution becomes multivalued, thus a weak solution with unclear physical relevance.

Therefore the existence of solitonic solutions only appears in a nonlinear and dispersive context, and is even in this context non-generic.

Related Physical Phenomena

We focus in this section on physical critical phenomena due to the nonlinearity in dispersive PDEs, mainly dispersive shock waves and blow-up.

Dispersive Shock Waves

We illustrate this notion by considering the KdV equation in which we introduce the slow variables x → x and t → t for a small parameter. We thus consider KdV in the following form

∂ t u + 6u∂ x u + 2 ∂ 3 x u = 0, x ∈ R, t ∈ R (1.26)
The limit → 0 is called the small dispersion limit. The corresponding dispersionless equation ( = 0) is the Hopf equation

∂ t + 6u∂ x u = 0. (1.27)
Its local classical solution built with the method of characteristics breaks down after a finite elapse of time. Thus a more general notion of solution has to be introduced.

Let u 0 ∈ L ∞ loc (R). A function u ∈ L ∞ loc (R × [0, +∞[) is called a weak solution to the Cauchy problem (1.27) if it satisfies +∞ -∞ +∞ 0 u ∂φ ∂t + u 2 2 ∂φ ∂x dtdx + +∞ -∞ u 0 (x)φ(x, 0)dx = 0 (1.28)
for all φ of class C 1 with compact support in R × [0, +∞[. One can then construct global weak, even entropy solutions3 to the Cauchy problem (1.27), in situations where no classical solution exists.

In fact, as soon as the gradient ∂ x u becomes important in the characteristic approach, the simplifications leading to a model as (1.27) are no longer valid. One has neglected terms in the original problems, which, even if they do not play a role where |∂ x u| is small, have an important impact where |∂ x u| is large. The zones where multivalued solutions appear correspond to the latter. Thus one has to use an enhanced model instead, and we describe below the two main possible approaches.

The first one consists in doing a dissipative regularization of the Hopf equation (1.27), naturally given by the viscous Burgers equation,

∂ t u + 6u∂ x u -ν∂ xx u = 0 (1.29)
with ν > 0. The term ν∂ xx u is called the dissipative term, and ν can represent, for example in fluid mechanics, the viscosity. Thanks to the entropy condition, the limiting solution to the viscous equation as the viscosity tends to zero is an entropy solution to (1.27).

With this model, solutions can be interpreted as discontinuity waves. We show the time evolution of the solution of (1. The second approach consists in considering a dispersive regularization to the Hopf equation (1.27), namely the KdV equation in the small dispersion regime (1.26).

With this model, the regularization of the breakup singularity occurs via the generation of rapid modulated oscillations in the region where the corresponding solution for Hopf has shocks, and the frequency of these oscillations increases as decreases. This phenomenon is called a dispersive shock, and a typical example is shown in Fig 1 .6 for the KdV equation, with initial data of the form u 0 = -sech 2 (x); see [START_REF] Grava | Numerical solution of the small dispersion limit of Korteweg de Vries and Whitham equations[END_REF] for more detail. An asymptotic description of these dispersive shock solutions to the KdV equation (1.26) can be given: the oscillatory behavior can be approximately described by the exact periodic solution of the KdV equation, where the spectral parameters are not constant and evolving according to the Whitham equations, see [START_REF] Whitham | Linear and nonlinear waves[END_REF] and [START_REF] Grava | Numerical solution of the small dispersion limit of Korteweg de Vries and Whitham equations[END_REF]. Gurevitch and Pitaevskii first proposed a mathematically consistent formulation of the problem of formation and evolution of a dispersive shock wave entirely in the framework of hyperbolic equations, see [START_REF] Gurevich | Nonstationary Structure of a collisionless Shock Wave[END_REF]:

Let the upper half-plane (t, x) be split into three domains:

{(t, x); t ≥ 0, x ∈ R} = {t ≥ 0; (-∞, x -(t)) ∪ [x -(t), x + (t)] ∪ (x + (t), ∞)}
Then, the asymptotic solution to equation (1.26) is governed by different equations.

For t < t c , it is described by the solution to the Hopf equation, and for t > t c , the oscillatory zone (roughly where the Hopf solution via the method of characteristics is multi-valued) is described by the solution of the Hopf equation for x / ∈ [x -(t), x + (t)]; and for x ∈ [x -(t), x + (t)] by the exact elliptic solution to KdV, with branch points dependent on (t, x) according to the Whitham equations (see [START_REF] Grava | Numerical solution of the small dispersion limit of Korteweg de Vries and Whitham equations[END_REF]).

The Lax-Levermore-Venakides construction, see [START_REF] Lax | The small Dispersion Limit of the Korteweg de Vries Equation I, II, III[END_REF], and [START_REF] Venakides | The Korteweg-de Vries Equation with small Dispersion: higher order Lax-Levermore Theory[END_REF], [START_REF] Venakides | The zero dispersion limit of the Korteweg de Vries equation for initial potential with nontrivial reflection coefficient[END_REF] can be considered as a formal justification of the validity of the direct Gurevich-Pitaevskii formulation of the problem in terms of the Whitham equations. The behavior of the solution of the KdV equation in the small dispersion limit, and the corresponding asymptotic solution can be seen in Fig. 1.7 for = 0.01. KdV and asymptotic solution with u 0 = -sech 2 (x) and = 0.01, see [START_REF] Grava | Numerical solution of the small dispersion limit of Korteweg de Vries and Whitham equations[END_REF].

The asymptotic description of these dispersive shocks for equations of NLS type is known only for particular cases. In this context, since the parameter has a similar role as the Planck constant in the Schrödinger equation, the limit → 0 is also called the semiclassical limit. For the (rescaled) one dimensional cubic NLS,

i ∂ t u + 2 ∂ xx u -2λ |u| 2 u = 0, (1.30)
Jin, Levermore and McLaughlin, see [START_REF] Jin | The behavior of solutions of the NLS equation in the semiclassical limit[END_REF], obtained the limit for the defocusing case by using inverse scattering theory. It is similar to the KdV case, dispersionless and Whitham equations are hyperbolic. A typical example is shown in Fig. 1.8, with u(0, x) = e -x 2 , = 0.5 and t max = 1. The initial pulse is broadened but gets steeper on both sides, before reaching the point of gradient catastrophe, where small oscillations appear. Only special initial data have been studied for the focusing cubic NLS equation in the semiclassical limit. Kamvissis, McLaughlin and Miller (see [START_REF] Kamvissis | Semiclassical Soliton Ensembles for the Focusing Nonlinear Schrödinger Equation[END_REF]) obtained the limit for real initial data, and Tovbis, Venakides and Zhou (see [START_REF] Tovbis | On semiclassical (zero Dispersion Limit) Solutions of the focusing nonlinear Schrödinger Equation[END_REF]) for initial data of the form u(0, x) = -sech x e (-iµ ln cosh x) . The dispersionless and the Whitham equations are elliptic in this case. Thus the solutions show a completely different behavior than in the defocusing case. A typical example is shown in Fig. 1.9, with u(0, x) = e -x 2 , = 0.1 and t max = 0.8. The square of the absolute value of u for the initial pulse is focused until it reaches its maximal height at time t = 0.26, where the strongest gradient appears. After this time the plot shows several smaller humps of similar shape as the one at breakup of the dispersionless equation.

Solutions to nonlinear dispersive PDEs can also have blow-up, i.e., a loss of regularity of the solution with respect to the initial data, as we illustrate in the next section.

Blow-up

The meaning of blow-up is context-dependent. In the differential equations setting, it usually refers to the fact that the solution of an evolution equation is defined only on a finite time domain. At the end of this interval, something 'bad' happens: either the solution goes to infinity, or it stops being smooth (in a way that makes the differential equation stop having sense). In all cases, that means that there is an explosion for a certain norm of the solution.

Figure 1.9: Solution to the focusing cubic NLS equation in semiclassical limit with u(0, x) = e -x 2 , = 0.1 and t max = 0.8.

We illustrate this notion with the very simple example:

∂ t u = -u + u 2 , u(t) ∈ R (1.31) Its solution is u(t) = 1 1-ce t . Hence, if c ≤ 0 (which corresponds to 0 < u(0) < 1)
, the solution exists globally in time for t ≥ 0, but if 0 < c < 1, (which corresponds to u(0) > 1), the solution blows up at t = log(1/c). We thus have global existence of solutions with small initial data and local existence of solutions with large initial data.

This behavior also occurs in many PDEs; for small initial data, linear damping terms can dominate the nonlinear terms, and global solutions exist. However for large initial data, the nonlinear blow-up overcomes the linear damping, and one only has local in time solutions. There are different ways in which blow-up can occur in PDEs, the solution itself may become unbounded (as in the ODEs setting), the explosion of a Sobolev norm for the solution, and/or a loss of regularity is observed, or the solution does not stay in a compact space, (for example lim u → +∞). This is blow-up in the narrow sense. On the other hand, spatial derivatives of the solution may become unbounded, due to the formation of some kind of singularity in the solution, as for example in equation (1.25) that we discussed in the previous section. This situation is also referred to a blow-up. In some cases, it is possible to continue a smooth solution past the blow-up time by a weak solution; in other cases, blow-up indicates a catastrophic breakdown in the ability of the PDE to model the original system.

For many of the PDEs under consideration, no conjecture exist to predict blow-up. The precise mechanism of blow-up remains an open problem, even if for some cases, one knows when blow-up can occur.

For example, consider the quintic NLS equation,

i ∂ t u + 2 ∆u -2λ |u| 4 u = 0.
(1.32)

In the focusing case (λ = -1), and with = 1, the solution to (1.32) blows up for any initial condition with negative energy, see [START_REF] Merle | On universality of blow-up profile for L 2 critical nonlinear Schrödinger equation[END_REF]. In the semiclassical regime ( → 0), we observe a blow-up of the solution, for any u 0 = 0, for sufficiently small. However, no dispersive shock waves appear before blow-up, as we can see in Fig. 1.10, where we show the solution of (1.32) for initial data u(0, x) = sech (x) and = 0.1. Another well known example for blow-up is the dimension critical generalized KdV equation

∂ t u + ∂ x ∂ xx u + u 5 = 0, (x, t) ∈ R × [0, T ) u(0, x) = u 0 (x), x ∈ R (1.33)
A complete review can be found for example in [START_REF] Merle | Existence of Blow-up Solutions in the Energy Space for the critical Generalized KdV Equation[END_REF] and for the generalized KP equations, see [START_REF] Klein | Numerical study of blow up and stability of solutions of generalized Kadomtsev-Petviashvili equations[END_REF].

Nonlinear effects in dispersive PDEs as dispersive shock waves and blow-up are mainly known for one-dimensional PDEs. No asymptotic description of dispersive shocks is known for (2 + 1)-dimensional PDEs, and no asymptotic description of the blow-up phenomena for the PDEs studied here. Efficient numerical simulations therefore become the only way to investigate these questions.

The KP and the DS II Equations, (2+1)-Dimensional PDEs

We focus in this work on the Kadomtsev-Petviashvili (KP) equations, and the Davey-Stewartson (DS) equations. Therefore we review some analytical aspects and explicit solutions in this section.

Analytic Properties of the KP Equations

The Lax representation of the KP equation (1.16) is given by

∂ xx Ψ + ∂ y Ψ = -(u + ρ)Ψ (1.34) ∂ t Ψ + 4∂ 3 x + 6u∂ x Ψ + 3∂ x uΨ = 3λ ∂ -1 x ∂ y u Ψ + αΨ (1.35)
where the solution u(x, y, t) of the KP equation plays the role of a scattering potential, Ψ(x, y, t, ρ) being the corresponding eigenfunction, ρ the spectral parameter, and α an arbitrary constant. The anti-derivative operator ∂ -1

x is defined as

∂ -1 x (f )(x) = 1 2 x -∞ f (x )dx - ∞ x f (x )dx (1.36)
and used in this form for IST 4 . This definition allows to write the KP equation in evolutionary form (i.e., the compatibility condition of (1.34) and (1.35) for u):

∂ t u + 6u∂ x u + ∂ xxx u + 3λ∂ -1 x ∂ yy u = 0, λ = ±1, (1.39)
which is equivalent, under suitable conditions of convergence and regularity, (including for example boundary conditions as periodic or rapidly decreasing at infinity), to the following system

∂ t u + 6u∂ x u + 2 ∂ xxx u + 3λ∂ y v = 0 u | t=0 = u I (x, y) ∂ y u = ∂ x v (1.40)
The validity of (1.39) requires that

+∞ -∞ ∂ yy u(x, y, t)dx = 0. (1.41)
4 It can also be viewed as the Fourier multiplier -i kx . In fact, since

± k 2 x + k 2 y = ±kx 1 + ky kx 2 (1.37) ≈ for ky kx 1 kx + 1 2 k -1 x k 2 y (1.38) we can infer that ± √ ∆ = ± ∂ 2 x + ∂ 2 y ≈ ∂x + 1 2 ∂ -1 x ∂yy
This condition imposes an infinite number of constraints on the initial datum, and even if this constraint is not satisfied at t = 0, it was shown in [START_REF] Fokas | The inverse spectral method for the KP I equation without the zero mass constraint[END_REF][START_REF] Molinet | Remarks on the mass constraint for KP type equations[END_REF] that the solution to the Cauchy problem associated to (1.39) will satisfy it for all t > 0 and that a discontinuity develops at t = 0; the resulting solution being thus not smooth in time. Numerical experiments in [START_REF] Klein | Numerical Study of oscillatory Regimes in the Kadomtsev-Petviashvili Equation[END_REF] indicate that the solution after an arbitrary small time step will develop an infinite 'trench' the integral over which just ensures that (1.41) is fulfilled.

The IST schemes for the KP equations are dependent on the value of λ = ±1, and they differ substantially from each other.

For the KP I equation (λ = -1), it was derived by Manakov [START_REF] Manakov | The inverse scattering transform for the time-dependant Schrödinger equation and Kadomtsev-Petviashvili equation[END_REF], Fokas and Ablowitz [START_REF] Fokas | On the inverse scattering of the time-dependent Schrödinger equation and the associated Kadomtsev-Petviashvili equation[END_REF], and Boiti, Leon and Pempinelli [START_REF] Boiti | Spectral transform and orthogonality relations for the Kadomtsev-Petviashvili I equation[END_REF], and is formulated in terms of a nonlocal Riemann-Hilbert problem. Unlike for KP I, the IST scheme for the KP II equation (λ = 1) cannot be formulated as a Riemann-Hilbert boundary value problem, since it turns out that the eigenfunctions are nowhere analytic. In this case, a generalization of the Riemann-Hilbert problem is required, the ∂ problem, first introduced by Beals and Coifman (see [START_REF] Beals | Scattering, transformations spectrales et équations d'évolution non linéaires[END_REF], [START_REF]Scattering, transformations spectrales et équations d'évolution non linéaire[END_REF]); and then used by Ablowitz, BarYaacov and Fokas (see [START_REF] Ablowitz | On the inverse scattering transform for the Kadomtsev-Petviashvili equation[END_REF]) to show that this approach is essential in the development of the IST scheme for KP II 5 .

For a recent review of the integrable aspects of KP see [START_REF] Fokas | Lax pairs: a novel type of separability[END_REF]. The complete integrability of the KP equations implies the existence of exact solutions. These include line solitons, i.e., solutions localized in one spatial direction and infinitely extended in another, for instance the 1-soliton given by

u(x, y, t) = 1 2 a 2 sech 2 1 2 a x -by - ωt a -x 0 , (a, b, x 0 ) ∈ R 3 , (1.42)
where ω is dependent on a, b, and given by the dispersion relation of the KP equation. These solutions are all unstable for large a for KP I, and for the KP II equation they are believed to be stable, although no formal proof exists.

The KPI equation has a two-dimensional solitary wave called a lump soliton. It is localized in all spatial directions, and given by

u(x, y, t) = 4 -(x + ay + 3(a 2 -b 2 )t) 2 + b 2 (y + 6at) 2 + 1/b 2 ((x + ay + 3(a 2 -b 2 )t) 2 + b 2 (y + 6at) 2 + 1/b 2 ) 2 , (a, b) ∈ R 2 (1.43)
It was found in 1977 by Manakov et al. [START_REF] Manakov | Two-Dimensional Solitons of the Kadomtsev-Petviashvili Equation and their Interaction[END_REF], and afterwards, various more general rational solutions to the KP I equation were obtained. Zaitsev [START_REF] Zaitsev | Formation of stationary waves by superposition of solitons[END_REF] found a solution which is localized in one direction and periodic in the second (a transformation of the form x → ix, y → iy exchanges these two directions). It has the form

u(ξ, y) = 2α 2 1 -β cosh(αξ) cos(δy) (cosh(αξ) -β cos(δy)) 2 , (α, β) ∈ R 2 , (1.44)
where

ξ = x -ct, c = α 2 4 -β 2 1 -β 2 , and δ = 3 1 -β 2 α 2 .
This solution is shown in Fig. 1.11, it is localized in x, periodic in y, and unstable as discussed in [START_REF] Klein | Numerical study of blow up and stability of solutions of generalized Kadomtsev-Petviashvili equations[END_REF]. No such real non-singular rational solution is known for KP II, but algebro-geometric solutions to the KP equations can be constructed on an arbitrary compact Riemann surface, see e.g. [START_REF] Dubrovin | Theta-functions and nonlinear equations[END_REF], [START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods: KdV and KP solutions[END_REF]. These solutions are in general almost periodic. Solutions on genus 2 surfaces, which are all hyperelliptic, are exactly periodic, but in general not in both x and y. A doubly periodic solution of KP II of genus 2 can be written as

u(x, y, t) = 2 ∂ 2 ∂x 2 ln θ (ϕ 1 , ϕ 2 ; B) (1.45)
where θ (ϕ 1 , ϕ 2 ; B) is defined by the double Fourier series

θ (ϕ 1 , ϕ 2 ; B) := ∞ m 1 =-∞ ∞ m 2 =-∞ e 1 2 m T Bm+im T ϕ (1.46)
where m T = (m 1 , m 2 ), and where B is a 2 × 2 symmetric, negative-definite Riemann matrix

B = b bλ bλ bλ 2 + d
, with real parameters λ = 0, b and d.

The phase variable ϕ has the form ϕ j = µ j x + ν j y + ω j t + ϕ j,0 , j = 1, 2. The solution can be seen in Fig. 1.12. It travels as the Zaitsev solution with constant speed in x-direction.

Analytic Properties of the DS Equations

We study here only the DS II equations (α = 1 in equations (1.10), (1.11)) since the elliptic operator for Φ can be inverted by imposing simple boundary conditions. For a hyperbolic operator acting on Φ boundary conditions for wave equations have to be used.

The Cauchy problem for the defocusing DS II equation, with decaying initial conditions, has a global in time solution, which disperses away the initial data after a long time, see [START_REF] Fokas | On the solvability of the N-wave, the Davey-Stewartson and the Kadomtsev-Petviashmli equation[END_REF]. The focusing DS II equation has a global solution under the condition that a certain norm of the initial data is small, see [START_REF] Fokas | On the solvability of the N-wave, the Davey-Stewartson and the Kadomtsev-Petviashmli equation[END_REF].

Theorem 1.1 If u 0 (x, y) ∈ S(R 2 )
, and if u 0 1 and u 0 ∞ are sufficiently small, then DS II has an unique solution u

∈ C 0 R, S(R 2 )
Without this small norm condition, no global existence result is known.

In fact it is known that DS II solutions can have blow-up: results by Sung [START_REF] Sung | Long-time decay of the solutions of the Davey-Stewartson II equations[END_REF] establish global existence in time for initial data u 0 ∈ L p , 1 ≤ p < 2 with a Fourier transform

F[u 0 ] ∈ L 1 ∩ L ∞ subject to the smallness condition ||F[u 0 ]|| L 1 ||F[u 0 ]|| L∞ < π 3 2 √ 5 -1 2 2 (1.47)
in the focusing case. There is no such condition in the defocusing case.

In [START_REF] Ozawa | Exact Blow-up Solutions to the Cauchy Problem for the Davey-Stewartson Systems[END_REF], Ozawa gave an example of an L 2 solution which blows up in finite time for the DS II system in the form

i∂ t u -∂ xx u + ∂ yy u = λ|u| 2 u + µφu ∂ xx φ + α∂ yy φ = |u| 2 xx (1.48) Theorem 1.2 (Ozawa solution) Let α > 0, µ = -2λ = 16α, ab < 0 and T = -a/b.
Denoting by u α (x, y, t) the function defined by

u α (x, y, t) = exp -i b 4(a + bt) (x 2 -y 2 ) v α (X, Y ) a + bt , (1.49) where v α (X, Y ) = 1 1 + α(X 2 + Y 2 ) , X = x a + bt and Y = y a + bt . (1.50)
Then, u α is a solution of (1.48) with

u α (t) 2 = v α 2 = π α (1.51) and |u α (t)| 2 → π α δ in S when t → T. (1.52)
where δ is the Dirac measure.

Beals and Coifman, see [START_REF]Multidimensional inverse scatterings and nonlinear partial differential equations[END_REF], and Wickerhauser, see [START_REF] Wickerhauser | Inverse scattering for the heat operator and evolutions in 2 + 1 variables[END_REF], gave a rigorous analytic treatment for scattering operators associated to DS II. Its complete integrability (see also [START_REF] Ablowitz | Nonlinear evolution equations in two and three dimensions[END_REF]) implies again that explicit solutions are known. Arkadiev and al, see [START_REF] Arkadiev | Inverse Scattering Transform Method and Soliton Solutions for Davey-Stewartson II Equation[END_REF], gave certain localized exact solutions to the DS II equation, known as lumps, by inverse scattering methods.

For example, the one-lump solution of DS II is written as

u(x, y, t) = 2 c1 (0) exp 2i(ξx -ηy + (ξ 2 -η 2 )t (x + 2ξt + z 1R (0)) 2 + (y + 2ηt + z 1I (0)) 2 + |c 1 (0)| 2 (1.53)
where c 1 and z 1 are parameters. It describes a localized nonsingular object which moves in the (x, y)-plane with the constant velocity (-2ξ, -2η) and decays as (

x 2 + y 2 ) -1 .
First observed by Gadylshin and Kyselev (see [START_REF] Gadyl Shin | Structural Instability of a Soliton for the Davey-Stewartson II Equation[END_REF]), the structural instability of one-lump solutions of the DS II equation has been proved by Pelinovski and Sulem, see [START_REF] Pelinovsky | Spectral decomposition for the Dirac System associated to the DSII Equation[END_REF]. They showed that for small data, localized initial conditions decay into radiation unless these data happen to be initial conditions corresponding exactly to lumps. That means that no other initial condition will eventually evolve into lumps, contrary to the KP case, see [START_REF] Klein | Numerical study of blow up and stability of solutions of generalized Kadomtsev-Petviashvili equations[END_REF].

Conclusion

We presented in this chapter the notion of dispersive PDEs, and properties of these PDEs, mainly complete integrability and critical physical phenomena under consideration in this work as dispersive shock waves and blow-up. We focused on dispersive PDEs in (2 + 1)-dimensions, the KP equation, and the DS II equation for which we want to study dispersive effects.

In view of the importance of these equations in applications and the open mathematical questions, efficient numerical algorithms are needed to enable extensive numerical studies. We describe the numerical treatment applied on these problems in the next chapter.

Chapter 2

Numerical Methods

The numerical solution of nonlinear dispersive PDEs is a highly complicated and problemdependent task due to the fact that the solution may contain dynamic interactions between shock waves, rarefaction waves and contact discontinuities. This implies that a method developed for a particular test problem may or may not work for another with stronger or weaker shocks and contact discontinuities; and that methods which are efficient in one space dimension may or not be easily extended to two or three space dimensions. The discussion in this chapter is about solving systems of conservation laws which are modeled by hyperbolic dispersive PDEs in one and higher (at least two) space dimensions. The numerical treatment is described, and based on the method of lines (MOL) see [START_REF] Hyman | The method of lines solution of partial differential equations[END_REF].

In the latter, we first select a method to discretize the differential equation in space and incorporate boundary conditions before integrating equations in time by using an appropriate method chosen in accordance with the properties of the spectrum of the discrete operator resulting from the discretization in space.

Space Discretization: Spectral Method

To describe numerically dispersive critical phenomena as discussed in Chap. 1, the use of efficient methods without introducing artificial numerical dissipation is necessary, in particular in the study of PDEs in several space dimensions. Typical numerical techniques used for the space discretization are finite elements, finite differences, and more recently, since forty years, spectral methods. Since the critical phenomena studied here are generally believed to be independent of the chosen boundary conditions, we study a periodic setting for simplicity, and use a spectral method based on Fourier series.

A Fourier Spectral Method

We describe in this part the Fourier spectral method that we use for the semi-discretization.

For the ease of the presentation, we restrict ourselves to the one space dimension. Spectral methods are a class of discretizations for differential equations which provide a way of translating an equation expressed in continuous space and time into a discrete equation which can be solved numerically, by approximating functions appearing in the differential equation via a sum of global, smooth and orthogonal functions. For a comprehensive review of spectral methods, see [START_REF] Fornberg | A practical Guide to pseudospectral Methods[END_REF], in which practical applications of collocation methods are presented and illustrated with examples and key Fortran code segments; see also [START_REF] Trefethen | Software, Environments, and Tools[END_REF], in which spectral collocation methods and their programming in Matlab are introduced and illustrated, and for a more theoretical review, see [START_REF] Canuto | Spectral methods, Scientific Computation[END_REF][START_REF]Evolution to complex geometries and applications to fluid dynamics[END_REF].

From the continuous to the discrete case.

The Fourier transform

F[u](k) of a function u(x) ∈ L 2 (R) is defined by F[u](k) = ∞ -∞ u(x)e -ikx dx (2.1)
where k ∈ R is the wavenumber. Conversely, u(x) can be reconstructed from F[u](k) by the inverse Fourier transform

u(x) = 1 2π ∞ -∞ F[u](k)e ikx dk, x ∈ R. (2.2)
Theorem 2.1 below makes a link between the smoothness of u and the corresponding decay of F[u].

Theorem 2.1 (see [START_REF] Trefethen | Software, Environments, and Tools[END_REF]):

Let u ∈ L 2 (R) have Fourier transform F[u].
1. If u has p -1 continuous derivatives in L 2 (R) for some p ≥ 0 and a p th derivative of bounded variation, then

F[u] = O(|k| -p-1 ) as |k| → ∞ (2.3)
2. If u has infinitely many continuous derivatives in L 2 (R), then

F[u] = O(|k| -m ) as |k| → ∞ (2.4)
for every m ≥ 0. The converse also holds.

3. If there exist α, c > 0 such that u can be extended to an analytic function in the complex strip |Imz| < α with u(• +iy) ≤ c uniformly for all y ∈ (-α, α), where

u(• +iy) ≤ c is the L 2 norm along the horizontal line Imz = y, then u α ∈ L 2 (R), where u α (k) = e α|k| F[u](k).
The converse also holds.

If u can be extended to an entire function and there exists

α > 0 such that |u(z)| = o e α|z| as |z| → ∞ for all complex values z ∈ C, then F[u] has compact support in [-α, α], that is, F[u](k) = 0, ∀ |k| > α. (2.5)
The converse also holds.

First Consequence of the Discretization

We consider a function u, 2π-periodic in x on the interval I = [0, 2π], and a regular subdivision {x 0 , x 1 , • • • , x N } of I with a step length h = 2π N and N even1 . The Discrete Fourier Transform (DFT) of u is defined by

vk = h N j=1 e -ikx j v j , k = - N 2 + 1... N 2 (2.6)
with v j = u(x j ), x j = 2jπ/N, for j = 1..N ; and the inverse discrete Fourier transform is given by

v j = 1 2π N 2 k=-N 2 +1
e ikx j vk , j = 1...N.

(2.7)

Due to the discretization of the spatial domain, the wavenumber domain is now the bounded interval

[ -π h , π h ] = [- N 2 , N 2 ]
This restriction is related to the phenomenon of aliasing, which refers to two kind of errors that can occur when a signal is sampled: different sampled signals become indistinguishable, or the signal reconstructed from samples differs from the original one. The error between the interpolating polynomial and the truncated series approximating u is called aliasing error.

The discrete Fourier transform (2.6) can also be expressed in terms of the Fourier transform of u:

Theorem 2.2 (Aliasing formula) Let u ∈ L 2 (R)
have a first derivative of bounded variation, and let v be the grid function on I defined by

v j = u(x j ). Then for all k ∈ [ -π h , π h ] vk = ∞ j=-∞ F[u] k + 2πj h (2.8)
Thus, the k-th mode of the interpolant of u depends on all modes of u, and the (k + jN )th wavenumber aliases the k-th wavenumber on the discrete grid of size N , they become indistinguishable.

As an illustration, we show the sampled functions of f 1 (x) = 2 cos(2πx + π/3) and f 2 (x) = 2 cos(2π6x + π/3) which are identical on the grid in Fig. 2.1. However, for well-resolved approximations, aliasing is not a real problem since the truncation and interpolation errors decay at the same rate, for sufficiently smooth function u; we have the discrete version of the theorem 2.1: Theorem 2.3 (see [START_REF] Trefethen | Software, Environments, and Tools[END_REF]): Let u ∈ L 2 (R) have a first derivative of bounded variation, and let v be the grid function on I defined by

v j = u(x j ). Then for all k ∈ [ -π h , π h ] 1. If u has p -1 continuous derivatives in L 2 (R)
for some p ≥ 0 and a p th derivative of bounded variation, then

|v(k) -F[u](k)| = O(h p+1 ) as h → 0 (2.9)

If u has infinitely many continuous derivatives in

L 2 (R), then |v(k) -F[u](k)| = O(h m ) as h → 0 (2.10)
for every m ≥ 0.

3. If there exist α, c > 0 such that u can be extended to an analytic function in the complex strip |Imz| < α with u(• +iy) ≤ c uniformly for all y ∈ (-α, α), then

|v(k) -F[u](k)| = O e -π(α-)/h as h → 0 (2.11)
for every > 0.

If u can be extended to an entire function, and if there exists

α > 0 such that |u(z)| = o e α|z| as |z| → ∞ for all complex values z ∈ C, then, provided h ≤ π α v(k) = F[u](k) (2.12)

Differentiation

In Fourier space, differentiation consists in multiplying each Fourier coefficient by ik,

where k is the corresponding wavenumber. More precisely, we have

u(x) = +∞ k=-∞ F[u] k e ikx , ( 2.13)
and for m > 0, the m th -derivative of u, is given by

∂ m x u(x) = +∞ k=-∞ (ik) m F[u] k e ikx ∈ L 2 (R) (2.14)
In physical space, differentiation is based on the values of the function u at the grid points

x j = 2πj N j=1..N .
We compute the values of the interpolation derivative, u j at the grid points (x j ) in terms of the values u l of the function at the same points, in the following way:

Evaluate the Discrete Fourier transform according to (2.6),

Perform the differentiation in Fourier space, as explained above, and then Transform back to physical space according the inversion formula (2.7).

This gives

(u ) j = 1 2π N 2 k=-N 2 +1
ŵk e ikx j , j = 1..N (2.15)

where ŵk = ikv k , for k = -N/2 + 1..N/2 -1, and ŵN/2 = 0 (2.16)
We compute in the same way the m th derivative of u j , by replacing (ik) by (ik) m in (2. [START_REF] Besse | Numerical Study of the Davey-Stewartson System[END_REF]), for all m > 0.

The analog of theorem 2.3 holds for u ∈ L2 (R) having a ν th derivative (ν ≥ 1) of bounded variation. Let w being the ν th spectral derivative of u on I, the estimates of theorem 2.3 hold uniformly for all x ∈ I, by replacing v by w j , F[u] by u (ν) (x j ) respectively, and O(p -1) by O(h p-ν ) in the first statement of the theorem.

A phenomenon closely related to the asymptotic behavior of the Fourier coefficients of a function that is related to the smoothness of that function 2 is known as the Gibbs' phenomenon. One can observe a particular behavior when a discontinuous function is expanded or interpolated with smooth functions, in the neighborhood of a point of discontinuity, as the appearence of oscillations. For example, consider a continuous function, and a discontinuous one,

f 1 : x → x if x ∈ [0, 1/2[ 1 -x if x ∈ [1/2, 1[ f 2 : x → 1 if x ∈ [0, 1/2[ -1 if x ∈ [1/2, 1[
The corresponding Fourier series are

S f 1 = 1 4 + n≥0 -cos(2π(2n + 1)x) π 2 (2n + 1) 2 and S f 2 = n≥0 4 sin(2π(2n + 1)x)) π(2n + 1)
The first is absolutely convergent, contrary to the second, and we show the corresponding behavior of the N -truncated Fourier series, which exhibits a Gibbs phenomenon for f 2 in Fig. 2.2.

Figure 2.2: The precision of the approximation depends on the smoothness of the function.

In practice, a higher resolution and a larger domain of computation can be used to avoid the difficulties associated with the Gibbs phenomenon.

Advantages of Spectral Methods with respect to alternative Approaches

To conclude this section, we list below the main benefits of using spectral methods to solve time dependent PDEs, with respect to Finite Differences (FD) or Finite Elements (FE) methods.

The first noticeable fact is that spectral discretizations, based for example on Fourier or Chebyshev series, provide very low error approximations. For analytic functions, these approaches can be exponentially convergent.

Secondly, spectral methods may require less memory and less time than FD or FE methods (in particular in several space dimensions). Indeed, the number of grid points necessary to obtain satisfactory accuracy can be very low and thus allows a time and memory efficient calculation. This 'memory-minimizing' property can be crucial in the implementation of parallel codes, because the memory used by an algorithm has a first-order effect on its runtime (time of completion) on modern parallel systems.

With the algorithm of Cooley and Tukey (1965) for the Fast Fourier transform (FFT) to compute the Discrete Fourier Transform, efficient algorithms for the transformation of bases for Fourier spectral methods exist and are implemented. Thus one does not need to implement these codes to develop a spectral method (most of the numerical libraries include a Fast Fourier Transform algorithm).

Spectral methods are nearly free of both dissipative or dispersive errors, which allows them to offer better accuracies and cost efficiencies in particular in simple geometries (boxes, spheres), but also for more complex shapes. They are of high performance in some areas (while the list is still growing) as for example turbulence modeling, weather prediction, nonlinear waves, seismic modeling, etc... (see [START_REF] Fornberg | A practical Guide to pseudospectral Methods[END_REF]). These approaches present also a surprising efficiency for many cases in which both solutions and variable coefficients are non-smooth or even discontinuous.

Thus, we proceed by approximating the spatial dependence via truncated Fourier series which leads for the studied equations to large systems of ODEs of the form

v t = Lv + N (v, t), ( 2.17)
where v denotes the (discrete) Fourier transform of u, and where L and N denote linear and nonlinear operators, respectively. The linear terms are uncoupled, and the diagonal matrix L is related to the wavenumbers. The resulting equations usually exhibit a stiffness property related to the linear part L, i.e., locally L contains widely separated eigenvalues. Their efficient solution requires the use of special integration schemes which allow relatively large time-step size without becoming unstable, as we will see in the next section.

Stiff Systems

In choosing the 'best' numerical method to integrate a system of ODEs, one has to consider accuracy, stability, storage requirements, computational complexity and the relative cost of the different methods.

For stiff equations of the form (2.17), the so-called stiffness property will lead to big restrictions on the choice of the time-integration scheme. We begin this section by introducing the notion of stiffness and related numerical difficulties for conventional methods.

Concept of Stiffness

There is no comprehensive and precise definition of the concept of stiffness. Therefore Lambert (see [START_REF] Lambert | Numerical Methods for ordinary differential Systems[END_REF]) considered stiffness as a phenomenon exhibited by a system, rather than a property of it. The latter is closely related to the notion of linear stability of a numerical method used to calculate the numerical solution of the following well-posed initial value problem

y = f (t, y), t ∈]0, T ], y(0) = a, ( 2.18)
where y ∈ C n , n being the dimension of the system.

To perform a stability analysis, we replace (2.18) by the problem

y = M (y -p(t)) + p (t), t ∈]0, T ] y(0) = a, ( 2.19)
where M = diag(m ii ), i = 1, .., n and m ii ∈ C. The general problem (2.18) can indeed be represented by equations of this form (2. [START_REF] Boiti | Spectral transform and orthogonality relations for the Kadomtsev-Petviashvili I equation[END_REF] with complex numbers M :

The behavior of solutions of (2.18) near a particular solution g(t) can be approximated by a Taylor series expansion which leads to

y = J(t, g(t))(y -g(t)) -f (t, g(t)) = J(t, g(t))(y -g(t)) -g (t), (2.20)
where J(t, g(t)) is the Jacobi matrix, which is assumed to be slowly varying in t, such that, locally, J(t, g(t)) can be taken as a constant matrix. For the ease of the exposition, we assume that J is diagonalizable 3 and we obtain the system (2. [START_REF] Boiti | Spectral transform and orthogonality relations for the Kadomtsev-Petviashvili I equation[END_REF]).

The solution of (2. [START_REF] Boiti | Spectral transform and orthogonality relations for the Kadomtsev-Petviashvili I equation[END_REF]) is y(t) = (a -p( 0))e M t + p(t).

(2.21)

The matrix M thus mainly influences the qualitative behavior of the solution. Three different cases are possible:

1. If Re(m ii ) 0 for all i then the problem is said to be unstable.

2. If Re(m ii ) > 0 and small for all i then the problem is said to be neutrally stable, and any conventional numerical method can be used.

3. If Re(m ii ) < 0 for all i and if there exists i and j, i = j such that Re(m jj ) Re(m ii ) is small, then the problem is stable, and the solution will tend to p(t) after a given time t called an initial transient.

Characterization of stiff problems:

A problem of the form (2.18) is said to be stiff if no solution component is unstable, or equivalently, if at least one eigenvalue of the Jacobian matrix has a negative part which is negative and large; and if the solution is slowly varying with respect to the negative real part of the eigenvalues of M . This characterization implies that a problem may be stiff in some intervals of t and not in others and as the term p(t) can suddenly change in dependence of time, stiff problems can exhibit various periods of rapid changes.

For a more complete review on the subject, we refer the reader to [START_REF] Hairer | Solving ordinary differential equations[END_REF] and [START_REF] Shampine | A user's view of solving stiff ordinary differential equations[END_REF].

Another way to characterize stiff problems found in the literature is the following:

Proposition 2.1 A numerical scheme is A-stable (or unconditionally A-stable) for all λ and ∆t if and only if

S ⊇ C - * := {z ∈ C, Re(z) < 0} (2.24)
Thus numerical methods used to integrate a stiff problem must have an A-stability region S extending to the whole left-half plane, thus avoiding stability restrictions on the time step, which can therefore be selected by accuracy conditions.

For the explicit Euler scheme, we have R(z) = (1 + z), which gives a A-stability region being the disk of center -1 and of radius 1. This scheme is thus A-stable if |1 + z| < 1, which imposes a restriction on the time step. We say that the scheme is conditionally A-stable.

For the implicit Euler scheme, R(z) = (1 -z) -1 , and so the A-stability region of this scheme is the set of z such that |1 -z| > 1. The scheme is thus unconditionally A-stable.

These examples are typical, explicit schemes can be only conditionally A-stable see [START_REF] Quarteroni | Numerical mathematics[END_REF], and implicit schemes are in general unconditionally A-stable.

Consider the following system,

x = M x, with M = -1 0 0 -µ (2.25)
with µ positive and large. Its solution is given by x(t) = x 1 (0)e -t , x 2 (0)e -µt , in which the second component tends to zero faster than the first one: There are two separated time scales in the system which is characteristic to stiff systems. If x 2 (0) is very small, or t large enough, one has x(t) ≈ (x 1 (0)e -t , 0). Using now an explicit Euler scheme to solve this problem,

x n+1 = (1 + ∆ t M )x n , we get x n = (x 1 (0)(1 -∆ t ) n , x 1 (0)(1 -µ∆ t ) n ).
Thus, to get a satisfying solution, ∆ t has to be chosen to satisfy the two following conditions

(1 -∆ t ) n ≈ e -n∆t (2.26) |1 -µ∆ t | < 1 (2.27)
The first condition is an accuracy condition, and the second a stability condition. Both impose a limitation on the time step length, and if µ is large, it is the stability condition ∆ t < 2 µ which limits the choice of the time step, instead of the aimed at accuracy of the numerical approximation. Conventional explicit numerical methods are thus inefficient in solving stiff systems.

There are some cases in which A-stability is not a desirable property, for example for methods with a rational stability function, for which the stiff components are damped out very slowly. Its strong requirement led to a restriction to the class of problems by the introduction of some other definitions of stability, less restrictive on the time step, as the A(α)-, the A 0 -and the L-stabilities, defined below.

Definition 2.3 A(α)-Stability

A method is said to be A(α)-stable, α ∈ (0, π/2), if

S ⊇ {z : -α < π -arg z < α} (2.28)
It is said to be A(0)-stable if it is A(α)-stable for some α ∈ (0, π/2).

Definition 2.4 A 0 -Stable A method is said to be A 0 -stable if S ⊇ {z : Re(z) < 0, Im(z) = 0} (2.29)
With these notions, Gear (1969) gave the following definition

Definition 2.5 Stiffly Stable

Let R 1 and R 2 be two regions of the complex plane defined as

R 1 := {z : Re(z) < -a}, R 2 := {z : -a < Re(z) < 0, -c < Im(z) < c},
where a and c are positive constants.

A method is said to be stiffly stable if

S ⊇ R 1 ∪ R 2 (2.
L-stability → A-stability → Stiff-stability → A(α)-stability → A(0)-stability → A 0 -stability
This hierarchy has implications for certain classes of methods that we discuss in the next section.

Choice of Time-Integration Schemes

In this section, we give results based on the stability hierarchy for the efficiency of Runge-Kutta and multistep methods in solving stiff systems, and select candidate methods for a numerical study of dispersive PDEs in (2 + 1)-dimensions.

Efficiency of Runge-Kutta and Multistep Methods in solving Stiff

Systems.

To construct high order methods from the Euler scheme, two different approaches are possible: the multistep methods, and the Runge-Kutta methods (RK). The main difference between them is that multistep methods are linear with respect to y n and f n = f (t n , y n ) and that they require only one evaluation of the function at each time step; their efficiency being thus increasable by increasing the number of time steps. In contrast, RK methods conserve the structure of one-step method, and their efficiency can be increased by multiplying the number of evaluations of the function at each time step, the linearity being thus sacrificed.

The general formulation of a s-stage Runge Kutta (RK) method for the initial value problem (2.18) is the following:

y n+1 = y n + ∆ t F (t n , y n , ∆ t ; f ), n ≥ 0
where F is an increment function defined by:

F (t n , y n , ∆ t ; f ) = s i=1 b i K i , ( 2.32
)

K i = f (t n + c i ∆ t , y n + ∆ t s j=1 a ij K j ), i = 1, 2...s, (2.33)
where s denotes the number of steps of the method, {(b i ) , (a ij ) , i, j = 1...s} are real numbers and

c i = s j=1 a ij .
The coefficients A = {a ij }, {c i } and {b i } characterize completely the method and can be represented in a Butcher table, see Table 2.1.

c 1 a 11 a 12 • • • a 1s c 2 . . . . . . • • • . . . . . . . . . . . . • • • . . . c s a s1 • • • • • • a ss b 1 b 2 • • • b s Table 2.

1: Butcher table of a RK scheme with s stages

If coefficients a ij are null for j ≥ i, i = 1..s, then each K i can be explicitly computed in dependence of the i -1 coefficients K 1 , ...., K i-1 already known. In this case, the method is explicit, otherwise it is implicit (denoted by IRK methods), and one has to solve a nonlinear system of dimension s to compute the K i . This increase of the computations makes the use of implicit schemes very expensive.

Key results for RK methods are recalled below (see [START_REF] Hairer | Solving ordinary differential equations[END_REF])

Theorem 2. [START_REF]On the numerical solution of the Sine-Gordon equation II. Performance of numerical schemes[END_REF] No explicit Runge-Kutta method is A-stable.

Proposition 2.3 (see [START_REF] Hairer | Solving ordinary differential equations[END_REF])

1. An IRK method has a polynomial stability function, R(z

) = P (z) Q(z) and is A-stable if and only if |R(iy)| ≤ 1 ∀y ∈ R (2.34)
and that R(z) is analytic for R(z) < 0. Table 2.2: Butcher tables of (from left to right): Gauss method (of order 4 with s = 2), Radau I and Radau II method (of order 3 with s = 2).

If the matrix

A = (a ij ) in
For multi-step methods, the result to point out is the following: Theorem 2.6 Dahlquist's second stability barrier (1963) An A-stable multi-step method cannot be more than second-order accurate.

This discussion of stiff systems and their efficient integration implies that implicit methods are recommended for stability reasons, but that means that nonlinear equations must be solved at each iteration, which represents an expensive computation compared to the obtained accuracy.

In addition, to avoid a pollution of the Fourier coefficients, and to allow the use of large time steps, we focus in the following on fourth-order schemes, and discuss some candidate methods for our purposes below.

Exponential Integrators

The construction and implementation of a class of methods called Exponential Integrators has been developed to study stiff and highly oscillatory semi-linear problems of the form

4 ∂ t u = Lu + N (u, t) (2.36)
where L represents the linear part responsible for the stiffness, and N the nonlinear part, which is slowly varying.

The first exponential integrators have been introduced by Certaine in 1960 (ETD2 and ETD3 based on Adams-Moulton methods) see [START_REF] Certaine | The Solution of ordinary differential Equations with large Time Constants[END_REF], followed by Lawson in 1967 (IF methods), see [START_REF] Lawson | Generalized Runge-Kutta Processes for stable Systems with large Lipschitz Constants[END_REF], as an alternative approach to deal with stiff systems. Since these pioneering works, the number of researchers working on this subject has steadily increased, leading to a large collection of exponential integrators. Amongst many others, one can distinguish the W-methods, and the Adaptive Runge-Kutta methods (ARK), see [START_REF] Strehmel | Behandlung steifer Anfangswertprobleme gewöhnlicher Differentialgleichungen mit adaptiven Runge-Kutta-Methoden[END_REF], [START_REF]B-Convergence Results for linearly implicit one Step Methods[END_REF], Integrating Factor methods (IF), [START_REF] Lawson | Generalized Runge-Kutta Processes for stable Systems with large Lipschitz Constants[END_REF], Exponential Time Differencing methods, see [START_REF] Nørsett | An A-stable Modification of the Adams-Bashforth Methods[END_REF] for ETD schemes based on Adams-Bashforth methods, and [START_REF] Friedli | Verallgemeinerte Runge-Kutta Verfahren zur Lösung steifer Differential gleichungssysteme[END_REF], for ETD based on explicit RK Methods.

Another class of methods known as Lie-group methods (LG) was first proposed in [START_REF] Munthe-Kaas | Lie-Butcher Theory for Runge-Kutta Methods[END_REF]. These transform the original differential equation to a new differential equation which evolves on a Lie algebra, and then take advantage of the fact that the Lie algebra is a linear space. Main examples of these methods are Crouch and Grossman (CG) methods, see [START_REF] Crouch | Numerical Integration of ordinary differential Equations on Manifolds[END_REF], Commutator-free Lie-Group (CF) methods, see [START_REF] Owren | Order Conditions for commutator-free Lie Group Methods[END_REF], and Runge-Kutta Munthe-Kaas (RKMK) methods, which also comprise the Generalized IF methods, see [START_REF] Krogstad | Generalized integrating Factor Methods for stiff PDEs[END_REF]. More recently, ETD schemes based on explicit RK methods up to order four were developed, see [START_REF] Cox | Exponential time differencing for stiff systems[END_REF], [START_REF] Krogstad | Generalized integrating Factor Methods for stiff PDEs[END_REF], and [START_REF] Hochbruck | Exponential Runge-Kutta Methods for semilinear parabolic Problems[END_REF].

This abundance of methods can be split into two categories, the exponential linear multistep methods and the exponential Runge-Kutta (multi-stage) methods; for a more comprehensive review on this topic, see [START_REF] Minchev | A review of exponential integrators for first order semi-linear problems[END_REF]. As we saw in the previous section that neither implicit methods nor multi-step methods are well suited for our problems, we only describe here explicit exponential Runge-Kutta integrators.

The main idea behind Exponential Integrators is to integrate exactly the linear part of the problem (which is mainly responsible for the stiffness), and then to use an appropriate approximation of the nonlinear part. Thus, the exponential function, and often functions closely related to it, appear in the formulation of the methods. Notice that the advances in computing approximations to the product of a matrix exponential with a vector played a big role in the renewed interest in the construction and implementation of exponential integrators.

The two main approaches to derive Exponential Integrator methods are the Integrating Factor Methods (IF), and the Exponential Time Differencing Methods (ETD).

Integrating Factor Methods

Lawson [START_REF] Lawson | Generalized Runge-Kutta Processes for stable Systems with large Lipschitz Constants[END_REF] suggested to improve the effect of the stiff linear part of equation (2.36) by using a change of the dependent variables (also called the Lawson transformation)

v(t) = e -Lt u(t) (2.37) Equation (2.36) becomes v (t) = e -Lt N (e Lt v(t), t) = e -Lt N (u(t), t), v(t = 0) = v 0 (2.38)
where v 0 = e -t 0 L u 0 . The designation "Integrating Factor methods" comes from the fact that if one applies directly the integrating factor e -tL to the equation ( 2.3. However, IFRK methods only work well on moderately stiff problems in which the solution tends to zero or is periodic, see [START_REF] Ehle | Generalized Runge-Kutta Processes for stiff Initial-Value Problems[END_REF].

Exponential Time Differencing Methods

The principle of Exponential Time Differencing Methods is the following: One uses the same integrating factor as in IF, (2.38), and integrates over a single time step of length ∆ t (i.e., between t n-1 and t n-1 + ∆ t ). One obtains the exact formula The approach now consists in approximating the nonlinear term N , and then solving the resulting integral.

u n+1 = e L∆t u n + e L∆t ∆t 0 e -Lτ N (u(t n + τ ), t n + τ ) dτ (2.40) c 1 c 2 a 21 e c 2 hL . . . . . . . . . c s a s1 e cshL • • • a ss-1 e (cs-c s-1 )hL b 1 e hL • • • b s-1 e (1-c s-1 )hL b s e (1-cs)hL
Explicit and implicit ETD schemes of arbitrary order have been first derived in [START_REF] Beylkin | A new Class of Time Discretization Schemes for the Solution of nonlinear PDEs[END_REF] in the following form:

Lemma : The exact solution of the initial value problem

u = Lu + N (u(t), t), u(t n-1 ) = u n-1 (2.41)
can be expressed in the form

u(t n-1 + ∆ t ) = e ∆tL u n-1 + ∞ k=0 ∆ k+1 t φ k+1 (∆ t L)N (k) n-1 (2.42)
where

N (k) n-1 = d i dt i t=t n-1
N (u(t), t) and where the φ k (z) are recursively defined as

φ 0 (z) = e z , φ k+1 (z) = φ k (z) -1 k! z , for k = 1, 2, ... (2.43)
However, this derivation is rather involved and does not give explicit formulas for the coefficients. Other derivations can be found in the literature, see [START_REF] Cox | Exponential time differencing for stiff systems[END_REF], [START_REF] Krogstad | Generalized integrating Factor Methods for stiff PDEs[END_REF] and references therein.

RKMK methods [START_REF] Munthe-Kaas | Lie-Butcher Theory for Runge-Kutta Methods[END_REF] are an example of ETD method based on classical fourth-order RK method, but are not well suited for problems where the stiffness is related to the linear term L in (2.36). The reason is that L is typically much larger than φ 1 (L) (see 2.43) and that the method involves the evaluation of products of the form LN (u).

High order ETDRK methods are difficult to derive and need to satisfy a large number of 'order conditions'. We present here three fourth-order ETDRK schemes which can be found in the literature: the fourth-order method of Cox and Matthews, see [START_REF] Cox | Exponential time differencing for stiff systems[END_REF], in Table 2.4, the fourth-order method of Krogstad, see [START_REF] Krogstad | Generalized integrating Factor Methods for stiff PDEs[END_REF], in Table 2.5, and the fourth-order method of Hochbruck and Ostermann (with five stages), see [START_REF] Hochbruck | Exponential Runge-Kutta Methods for semilinear parabolic Problems[END_REF], in Table 2.6.

0 1 2 1 2 φ 1,2 1 2 0 1 2 φ 1,3 1 1 2 φ 1,3 (φ 0,3 -1) 0 φ 1,3 φ 1 -3φ 2 + 4φ 3 2φ 2 -4φ 3 2φ 2 -4φ 3 -φ 2 + 4φ 3 Table 2.4: Fourth order method of Cox-Matthews 0 1 2 1 2 φ 1,2 1 2 1 2 φ 1,3 -φ 2,3 φ 2,3 1 φ 1,4 -2φ 2,4 0 2φ 2,4 φ 1 -3φ 2 + 4φ 3 2φ 2 -4φ 3 2φ 2 -4φ 3 -φ 2 + 4φ 3
Table 2.5: Fourth order method of Krogstad

In [START_REF] Hochbruck | Exponential Runge-Kutta Methods for semilinear parabolic Problems[END_REF], Hochbruck and Ostermann showed that in general, the construction of stiff fourth order exponential Runge-Kutta method requires at least five stages, without which, the method suffers from an order reduction, due to not satisfying some of the stiff order conditions derived in this paper. They showed that the fourth order Integrating Factor RK method (IRK4, see Table 2.2) has only stiff order one, that the fourth order method of Cox and Matthews (see Table 2.4 and [START_REF] Cox | Exponential time differencing for stiff systems[END_REF]) has only stiff order two and the fourth order method of Krogstad (see Table 2.5 and [START_REF] Krogstad | Generalized integrating Factor Methods for stiff PDEs[END_REF]) has a stiff order three in the worst case. For each of these methods, higher order convergence is possible, if additional smoothness conditions are satisfied by the problem .

Both Cox-Matthews' and Krogstad's schemes are, however, four-stage methods, whereas the Hochbruck-Ostermann method is a five-stage method that has stiff order four. Thus 0

1 2 1 2 φ 1,2 1 2 1 2 φ 1,3 -φ 2,3 φ 2,3 1 φ 1,4 -2φ 2,4 φ 2,4 φ 2,4 1 2 1 2 φ 1,5 -2a 5,2 -a 5,4 a 5,2 a 5,2 1 4 φ 2,5 -a 5,2 φ 1 -3φ 2 + 4φ 3 0 0 -φ 2 + 4φ 3 4φ 2 -8φ 3 with a 5,2 = 1 2 φ 2,5 -φ 3,4 + 1 4 φ 2,4 - 1 2 φ 3,5 , a 5,4 = 1 4 φ 2,3 -a 5,2 (2.44)
and φ i,j = φ i (-c j hL), j = 2..s.

(2.45) Table 2.6: Fourth order method of Hochbruck and Ostermann all these methods should show the same convergence rate in the non-stiff regime of the equation, but could differ for some problems in the stiff regime. We use in the following this characterization to detect these different regimes in the considered problem.

Notice that these results [START_REF] Hochbruck | Exponential Runge-Kutta Methods for semilinear parabolic Problems[END_REF] were established for stiff semi-linear parabolic problems, and that the applicability to hyperbolic PDEs of the type studied in this work is not obvious. One of the purposes of our studies is to get some experimental insight whether the Hochbruck-Ostermann theory holds also in this case. As the numerical tests show in the following, order reduction can be observed in some cases, but not in all, as it was the case in the experiments of Kassam and Trefethen, see [START_REF] Kassam | Fourth-Order Time-Stepping for stiff PDEs[END_REF].

The exponential methods described above are of multi-stage type, requiring only one previous evaluation of the nonlinear term N . Such methods are convenient to use, and also typically have the advantage of smaller error constants and larger stability regions than multi-step methods.

The main technical problem in the use of ETD schemes is the efficient and accurate numerical evaluation of the functions

φ i (z) = 1 (i -1)! 1 0 e (1-τ )z τ i-1 dτ, i = 1, 2, 3, 4,
i.e., functions of the form (e z -1)/z and higher order generalizations thereof, where one has to avoid cancellation errors. Kassam and Trefethen [START_REF] Kassam | Fourth-Order Time-Stepping for stiff PDEs[END_REF] used complex contour integrals to compute these functions.

The approach is straight forward for diagonal operators L that occur in our case because of the use of Fourier methods: one considers a unit circle around each point z and computes the contour integral with the trapezoidal rule which is known to be a spectral method in this case. Schmelzer [START_REF] Schmelzer | The fast evaluation of matrix functions for exponential integrators[END_REF] made this approach more efficient by using the complex contour approach only for values of z close to the pole, e.g. for |z| < 1/2. For the same values of z the functions φ i can be computed via a Taylor series. These two independent and very efficient approaches allow a control of the accuracy. We find that just 16 Fourier modes in the computation of the complex contour integral are sufficient to determine the functions φ i to the order of machine precision. Thus we avoid problems reported in [START_REF] Berland | Solving the nonlinear Schödinger Equation using exponential Integrators[END_REF], where machine precision could not be reached by ETD schemes due to inaccuracies in the determination of the φ i -functions. The computation of these functions takes only negligible time for the 2 + 1-dimensional equations studied in this work, especially since it has to be done only once during the time evolution. We find that ETD as implemented in this way has the same computational costs as the other used schemes.

Other possible Approaches

We finally discuss here several other approaches proposed in the literature to solve efficiently equations of the form (2.36) with a linear stiff part: IMEX methods in the form Driscoll's composite Runge-Kutta method, symplectic integrators, and time splitting methods.

IMEX Methods

The idea of IMEX methods (see e.g. [START_REF] Chan | Fourier methods with extended stability intervals for the Korteweg de Vries equation[END_REF] for KdV) is the use of a stable implicit method for the linear part of equation (2.17) and an explicit scheme for the nonlinear part which is assumed to be non-stiff. In [START_REF] Kassam | Fourth-Order Time-Stepping for stiff PDEs[END_REF] such schemes did not perform satisfactorily for dispersive PDEs which is why we only consider a more sophisticated variant here. Fornberg and Driscoll [START_REF] Driscoll | A fast spectral Algorithm for nonlinear Wave Equations with linear Dispersion[END_REF] provided an interesting generalization of IMEX by splitting also the linear part of the equation in Fourier space into regimes of high, medium, and low frequencies, and by using adapted numerical schemes in each of them. They considered the NLS equation as an example.

Driscoll's composite Runge-Kutta Method

Driscoll's [START_REF] Driscoll | A composite Runge-Kutta Method for the spectral Solution of semilinear PDEs[END_REF] idea was to split the linear part of the equation in Fourier space just into regimes of high and low frequencies. He used the explicit fourth order RK integrator for the low frequencies and the linearly implicit RK method of order three for the high frequencies. He showed that this method is in practice of fourth order over a wide range of step sizes. We confirm this in our experiments for the cases where the method converges, which it fails to do, however, sometimes for very stiff problems (see chapter 3). In particular, he used this method for the KP II equation for the two phase solution (1.12) we will also discuss in the next chapter. We call the method DCRK in the following.

Symplectic Integrators

Symplectic integrators are used for Hamiltonian systems of ODEs, in particular for systems of ODEs obtained after a semi-discretization of Hamiltonian PDEs by a symplecticconservative scheme (as pseudospectral method for example). The main problem with symplectic integrators is that, to the best of our knowledge, explicit scheme can be constructed if and only if the Hamiltonian system associated to the PDE is separable, and this is not the case for all dispersive PDEs. In other cases, symplectic schemes are known to be highly nonlinear and implicit, thus they are computationally expensive in comparison to traditional schemes, see [START_REF]On the numerical solution of the Sine-Gordon equation II. Performance of numerical schemes[END_REF], [START_REF] Islas | Geometric integrators for the nonlinear Schrödinger equation[END_REF]. On the other hand, even if symplectic integrators provide an 'expensive' but valuable tool for studies of long-time behavior of nonlinear PDEs, (and in particular they have been proven to be both accurate and efficient in the long-time approximation of solutions to Hamiltonian ODEs), both high-order spatial approximations and large values of resolution seem to reduce or eliminate these advantages, see [START_REF]On the numerical solution of the Sine-Gordon equation II. Performance of numerical schemes[END_REF], [START_REF] Ablowitz | On the numerical solution of the Sine-Gordon equation I. Integrable discretizations and homoclinic manifolds[END_REF] and [START_REF] Islas | Geometric integrators for the nonlinear Schrödinger equation[END_REF].

To confirm these facts, and since some implicit Runge-Kutta methods are symplectic and widely recommended for the case of a non-separable Hamiltonian, as the 2-Gauss method (IRK4 see Table 2.2), we will include IRK4 in some cases in our experiments.

More precisely, the implicit character of this method requires the iterative solution of a high dimensional system at every step which is done via a simplified Newton method. For the studied examples in the form (2.36), we have to solve equations of the form y = Ay + b(y) for y, where A is a linear operator independent of y, and where b is a vector with a nonlinear dependence on y. These are solved iteratively in the form

y n+1 = (1 -A) -1 b(y n ).
By treating the linear part that is responsible for the stiffness explicitly as in an IMEX scheme, the iteration converges in general quickly. Without taking explicit care of the linear part, convergence will be extremely slow. The iteration is stopped once the L ∞ norm of the difference between consecutive iterates is smaller than some threshold (in practice we work with a threshold of 10 -8 ). Per iteration the computational cost is essentially 2 FFT/IFFT pairs. Thus the IRK4 scheme can be competitive with the above explicit methods which take 3 or 4 FFT/IFFT pairs per time step if not more than 2-3 iterations are needed per time step. This can happen in the experiments of Chapter 3, in the nonstiff regime, but is not the case in the stiff regime. We only test this scheme where its inclusion appears interesting and where it is computationally not too expensive.

Splitting Methods

Splitting methods are convenient if an equation can be split into two or more equations which can be directly integrated. The motivation for these methods is the Trotter-Kato formula [START_REF] Trotter | On the Product of semi-Groups of Operators[END_REF] lim n→∞ e -tA/n e -tB/n n = e -t(A+B) (2.46) where A and B are certain unbounded linear operators, and either t ∈ iR or t ∈ R, t ≥ 0, and A and B are bounded from above. In particular this includes the cases studied by Bagrinovskii and Godunov in [START_REF] Bagrinovskii | Difference schemes for multi-dimensional problems[END_REF] and by Strang in [START_REF] Strang | On the Construction and Comparison of Difference Schemes[END_REF]. For hyperbolic equations, first references are Tappert [START_REF] Tappert | Numerical Solutions of the Korteweg-de Vries Equation and its Generalizations by the split-step Fourier Method[END_REF] and Hardin and Tappert [START_REF] Hardin | Applications of the split-step Fourier method to the numerical solution of nonlinear and variable coefficient wave equations[END_REF] who introduced the split step method for the NLS equation.

The idea of these methods for an equation of the form

∂ t u = (A + B) u
is to write the solution in the form

u(t) = exp(c 1 tA) exp(d 1 tB) exp(c 2 tA) exp(d 2 tB) • • • exp(c k tA) exp(d k tB)u(0)
where (c 1 , . . . , c k ) and (d 1 , . . . , d k ) are sets of real numbers that represent fractional time steps. Yoshida [START_REF] Yoshida | Construction of higher order symplectic integrators[END_REF] gave an approach which produces split step methods of any even order.

Conclusion

Our approach for the numerical treatment of dispersive PDEs is based on the Method of Lines, using a Fourier spectral method with periodic boundary conditions for the spatial discretization. This not only gives spectral accuracy in the spatial coordinates, but also minimizes the introduction of numerical dissipation which is important in the study of dispersive effects. The resulting systems of ODEs are classical examples of stiff equations where the stiffness is related to the linear part and which consequently, as discussed in section 2.2, require the use of special time-integration schemes. To solve our problems which deals with stiff systems and rapid oscillations, we selected some time-integration schemes as promising candidates. As already noticed a method working well for a particular stiff test problem may not work for another, which means that the most efficient method has to be identified for each PDE, but also for different regimes in which it is studied (mainly stiff or not stiff). This is the purpose of the next chapter, in which we compare the performance of the above selected numerical schemes (several explicit fourth order schemes mainly related to exponential integrators for various examples in a similar way as in the work by Kassam and Trefethen [START_REF] Kassam | Fourth-Order Time-Stepping for stiff PDEs[END_REF] and in [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger Equation[END_REF] for KdV and NLS) for the study in both stiff and non-stiff regimes of the Kadomtsev-Petviashvili equation, and of the Davey-Stewartson II equation.

Chapter 3

Comparison of Time Stepping Schemes

Introduction

Purely dispersive partial differential equations as the Korteweg-de Vries equation, the nonlinear Schrödinger equation and higher dimensional generalizations thereof can have solutions which develop a zone of rapid modulated oscillations in the region where the corresponding dispersionless equations have shocks or blow-up. Although an asymptotic description of these dispersive shocks is well known for certain integrable PDEs as KdV [START_REF] Lax | The small Dispersion Limit of the Korteweg de Vries Equation I, II, III[END_REF][START_REF] Venakides | The Korteweg-de Vries Equation with small Dispersion: higher order Lax-Levermore Theory[END_REF][START_REF] Deift | New result in small dispersion KdV by an extension of the steepest descent method for Riemann-Hilbert problems[END_REF] and the NLS equation for certain classes of initial data [START_REF] Jin | The behavior of solutions of the NLS equation in the semiclassical limit[END_REF][START_REF] Kamvissis | Semiclassical Soliton Ensembles for the Focusing Nonlinear Schrödinger Equation[END_REF][START_REF] Tovbis | On semiclassical (zero Dispersion Limit) Solutions of the focusing nonlinear Schrödinger Equation[END_REF], no such description is known for (2 + 1)-dimensional PDEs which, in addition, can have solutions which blow up (see Chap. 1). To numerically study such phenomena, fourth order timestepping in combination with spectral methods are beneficial to resolve the steep gradients in the oscillatory region (see Chap. 2). We compare in this chapter the performance of several fourth order time integration schemes discussed in the previous chapter for the KP and the DS II equation.

The KP equation is written here in the form

∂ x ∂ t u + 6u∂ x u + 2 ∂ xxx u + λ∂ yy u = 0, λ = ±1, ( 3.1)
where (x, y, t) ∈ R x × R y × R t and where 1 is a small scaling parameter introduced as in the KdV case 1 . The limit → 0 is the dispersionless limit. Higher dimensional generalizations of the KP equations, where the derivative ∂ yy is replaced by the Laplacian in the transverse coordinates, ∆ ⊥ = ∂ yy + ∂ zz , are important for instance in acoustics. The numerical problems to be expected there are the same as in the 2 + 1-dimensional case studied here.

The Davey-Stewartson II equation can be written in the form

i ∂ t u + 2 ∂ xx u -2 ∂ yy u + 2ρ Φ + |u| 2 u = 0, ∂ xx Φ + ∂ yy Φ + 2 |u| 2 xx = 0, (3.2)
where ρ takes the values ±1, where 1 is again a small dispersion parameter, and where Φ is a mean field. Since has the same role as the in the Schrödinger equation, the limit → 0 is also called the semiclassical limit in this context.

Both KP and DS II being completely integrable by IST, many explicit solutions are known and represent popular test cases for numerical algorithms. But as we will show in the example of KP, these exact solutions, typically solitons, often test the equation in a regime where stiffness is not important. The main challenge in the study of critical phenomena as dispersive shocks and blow-up is, however, the numerical resolution of strong gradients in the presence of which the above equations are stiff. As already noted this has the important consequence that algorithms that perform well for solitons might not be efficient in the stiff regime of the corresponding PDE.

In this chapter we are mainly interested in the numerical study of the KP and the DS II equations for Schwartzian initial data in the small dispersion limit. The latter implies that we can treat the problem as essentially periodic, and that we can use Fourier methods. After spatial discretization we thus face a system of ODEs of the form (2.17). Since we need to resolve high wavenumbers, these systems will be in general rather large. The PDEs studied here have high order derivatives in the linear part L of (2.17), whereas the nonlinear part N contains only first derivatives. This means that the stiffness in these systems is due to the linear part. The latter will thus be treated with adapted methods detailed in the previous chapter, whereas standard methods can be used for the nonlinear part. We restrict the analysis to moderate values of the dispersion parameter to be able to study the dependence of the different schemes on the time step in finite CPU time. For smaller values of see for instance [START_REF] Klein | Numerical Study of oscillatory Regimes in the Kadomtsev-Petviashvili Equation[END_REF] for the KP equation.

We compare several numerical schemes for equations of the form (2.17) as in [START_REF] Kassam | Fourth-Order Time-Stepping for stiff PDEs[END_REF] and [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger Equation[END_REF]. The PDEs are studied for a periodic setting with periods 2πL x and 2πL y in x and y respectively. We give the numerical error in dependence of the time step N t as well as the actual CPU time as measured by MATLAB (all computations are done on a machine with Intel 'Nehalem' processors with 2.93 GHz with codes in MATLAB 7.10). The goal is to provide some indication on the actual performance of the codes in practical applications. Since MATLAB is using in general a mixture of interpreted and precompiled embedded code, a comparison of computing times is not unproblematic. However, it can be done in the present context since the main computational cost is due to two-dimensional fast Fourier transformations (FFT).

For the KP equations all considered approaches (with the exception of the Hochbruck-Ostermann ETD scheme which uses 8 FFT commands per time step) use 6 (embedded) FFT commands per time step as was already pointed out in [START_REF] Driscoll | A composite Runge-Kutta Method for the spectral Solution of semilinear PDEs[END_REF]. Note that an additional FFT/IFFT pair is needed per time step since the method is set up in Fourier space, and since we are interested in the solution in physical space. For the DS II equation, the above numbers are doubled since the computation of the mean field Φ takes another FFT/IFFT pair per intermediate step. The φ-functions in the ETD schemes are also computed via FFT. It can be seen that this can be done with machine precision in a very efficient way. Since the φ-functions have to be obtained only once in the computation and since the studied problems are computationally demanding, this only has a negligible effect on the total CPU time in the experiments. The numerical error is the L 2 -norm of the difference of the numerical solution and an exact or reference solution, normalized by the L 2 -norm of the initial data. It is denoted by ∆ 2 .

KP Equations

In this section we study the efficiency of several numerical schemes in solving Cauchy problems for the KP equations. We first review previous numerical studies of the KP equations, and then select exact solutions to test numerical codes. We then compare the performance of the codes for the exact solutions and for a typical example in the small dispersion limit.

First numerical studies of exact solutions to the KP equations were performed in [START_REF] Wineberg | Implicit spectral Methods for Wave Propagation Problems[END_REF] and [START_REF] Xu | Local discontinuous Galerkin Methods for two Classes of two-dimensional nonlinear Wave Equations[END_REF]. A complete mathematical description of the small dispersion limit being not yet achieved (see [START_REF] Bogaevskiȋ | On Korteweg-de Vries, Kadomtsev-Petviashvili, and Boussinesq Equations in the Theory of Modulations[END_REF], [START_REF] Krichever | The averaging Method for two-dimensional "integrable[END_REF]), Klein, Markowich and Sparber [START_REF] Klein | Numerical Study of oscillatory Regimes in the Kadomtsev-Petviashvili Equation[END_REF] investigated the problem numerically. They provided numerical evidences to the fact that the corresponding dispersionless equation (dKP),

∂ x (∂ t u + u∂ x u) + λ∂ yy u = 0 λ = ±1 (3.3)
develops shocks in finite time (similarly to the KdV case), using a dissipative regularization of the dKP equation and that the solutions to the dKP equation give the correct limiting behavior of (3.1) as tends to 0 before breakup. Oscillatory zones are observed in the solution to (3.1) for small, but still non-zero , certainly due to the dispersive term ∂ xxx u in (3.1) which smooths out shocks of the solutions to (3.3), and the frequency of these oscillations increases as decreases, as in the KdV case.

In this section, we look for KP solutions that are periodic in x and y, i.e., for solutions on T 2 × R. This includes for numerical purposes the case of rapidly decreasing functions in the Schwartz space S(R 2 ) if the periods are chosen large enough that |u| is smaller than machine precision (we work with double precision throughout the chapter) at the boundaries of the computational domain. Notice, however, that solutions to Cauchy problems with Schwartzian initial data u 0 (x, y) will not stay in S(R 2 ) unless u 0 (x, y) satisfies an infinite number of constraints see (1.41). This behavior can be already seen on the level of the linearized KP equation, see e.g. [START_REF] Boiti | Solutions of the KPI equation with smooth initial data[END_REF][START_REF] Klein | Numerical Study of oscillatory Regimes in the Kadomtsev-Petviashvili Equation[END_REF], where the Green's function implies a slow algebraic decrease in y towards infinity. This leads to the formation of tails with an algebraic decrease to infinity for generic Schwartzian initial data. The amplitude of these effects grows with time (see for instance [START_REF] Klein | Numerical Study of oscillatory Regimes in the Kadomtsev-Petviashvili Equation[END_REF]). In our periodic setting this will give rise to echoes and a weak Gibbs phenomenon at the boundaries of the computational domain.

The latter implies that we cannot easily reach machine precision as in the KdV case unless we use considerably larger domains. As can be seen from computations in the small dispersion limit below and the Fourier coefficients in Sect. 3.4, we can nonetheless reach an accuracy of better than 10 -10 on the chosen domain. For higher precisions and larger values of t, the Gibbs phenomena due to the algebraic tails become important. To avoid the related problems, we always consider initial data that satisfy (1.41). A possible way to achieve this is to consider data that are x-derivatives of periodic or Schwartzian functions.

The most popular exact KP solutions are line solitons (1.42) discussed in Chap. 1, but such solutions typically have an angle not equal to 0 or 90 degrees with the boundaries of the computational domain, which leads to strong Gibbs phenomena. This implies that these solutions are not a good test case for a periodic setting. If the angle is 0 or 90 degrees, the solution only depends on one of the spatial variables and thus does not test a true 2d code. The lump soliton for KP I (1.43) is localized in all spatial directions, but only with algebraic fall off towards infinity. This would again lead to strong Gibbs phenomena in our setting. The Zaitsev solution (1.44) to KP I and the doubly periodic solution to KP II on an hyperelliptic Riemann surface of genus 2 (1.12) are thus preferred here due to their periodicity properties.

The 2-dimensional Fourier transform of u is written here in the form

F[u] := û(k x , k y , t) = R 2
u(x, y, t)e -ikxx-ikyy dxdy.

(3.4)

The KP equation is discretized first in space by Fourier spectral methods, and we get an equation of the form (2.17) with

L = -i λ k 2 y k x -2 k 3 x (3.5) and N (v, t) = -3ik x F F -1 (v) 2 (3.6)
where v denotes the 2-dimensional discrete Fourier transform of u.

We compare the efficiency of the following schemes: Cox and Matthews' ETD scheme (see Table 2.4 and [START_REF] Cox | Exponential time differencing for stiff systems[END_REF]), Krogstad's ETD scheme (see Table 2.5 and [START_REF] Krogstad | Generalized integrating Factor Methods for stiff PDEs[END_REF]), Hochbruck and Ostermann ETD scheme (see Table 2.6 and [START_REF] Hochbruck | Exponential Runge-Kutta Methods for semilinear parabolic Problems[END_REF]), the IFRK4 scheme (see Table 2.3) and Driscoll's composite Runge-Kutta scheme (see [START_REF] Driscoll | A composite Runge-Kutta Method for the spectral Solution of semilinear PDEs[END_REF]).

As indicated in the previous chapter, the use of IRK (see Table 2.2) can be competitive only in the non-stiff regime, and so is tested only in these cases.

The KP equation can be split into the following equations,

∂ t u + 6u∂ x u = 0 (3.7) (F[u]) t -ik 3 x F[u] + λ ik 2 y k x F[u] = 0 (3.8)
The Hopf equation (3.7) can be integrated in implicit form with the method of characteristics (see chapter 1), and the linear equation in Fourier space (3.8) can be directly integrated, but the implicit form of the solution of the former makes an iteration with interpolation to the characteristic coordinates necessary that is computationally too expensive. Therefore we will not use splitting methods here.

Remark 3.1

The standard 4th order Runge-Kutta scheme did not converge for any of the studied examples for the used time steps. The reason is that the Fourier multiplier -i/k x imposes very strong stability restrictions on the scheme.

Numerical Solution of Cauchy Problems for the KP I Equation

Zaitsev solution

We first study the case of the Zaitsev solution (1.44) with α = 1 and β = 0.5. Notice that this solution is unstable against small perturbations as shown numerically in [START_REF] Klein | Numerical study of blow up and stability of solutions of generalized Kadomtsev-Petviashvili equations[END_REF], but that it can be propagated with the expected numerical precision by the used codes. As initial data we take the solution centered at -L x /2 (we use L x = L y = 5) and propagate it until it reaches L x /2. The computation is carried out with 2 11 × 2 9 points for (x, y) ∈ [-5π, 5π] × [-5π, 5π] and t ≤ 1. The decrease of the numerical error is shown in Fig. 3.1 in dependence of the time step and in dependence of the CPU time. A linear regression analysis in a double logarithmic plot (log 10 ∆ 2 = -a log 10 N t + b) is presented in Fig. 3.1, where we can see that all schemes show a fourth order behavior. We find a = 4.32 for the Integrating Factor method, a = 4.38 for DCRK method, a = 3.93 for Krogstad's ETD scheme, a = 4 for the Cox-Matthews scheme, and a = 3.98 for the Hochbruck-Ostermann scheme.

In this context the DCRK method performs best, followed by the ETD schemes that have almost identical performance (though the Hochbruck-Ostermann method uses more internal stages and thus more CPU time in Fig. 3.1). It can also be seen that the various schemes do not show the phenomenon of order reduction as discussed in [START_REF] Hochbruck | Exponential Runge-Kutta Methods for semilinear parabolic Problems[END_REF], which implies that the Zaitsev solution tests the codes in a non-stiff regime of the KP I equation.

Small dispersion limit for KP I

To study KP solutions in the limit of small dispersion ( → 0), we consider Schwartzian initial data satisfying the constraint (1.41). As in [START_REF] Klein | Numerical Study of oscillatory Regimes in the Kadomtsev-Petviashvili Equation[END_REF] we use data of the form By numerically solving the dispersionless KP equation (3.3), we determine the critical time of the appearance of a gradient catastrophe by the breaking of the code, see [START_REF] Klein | Numerical Study of oscillatory Regimes in the Kadomtsev-Petviashvili Equation[END_REF].

u 0 (x, y) = -∂ x sech 2 (R) where R = x 2 + y 2 . ( 3
To study dispersive shocks, we run the KP codes for some time larger than this critical time. The solution can be seen in Fig. 3.2. It develops tails with algebraic fall off towards infinity. The wave fronts steepen on both sides of the origin. In the regions of strong gradients, rapid modulated oscillations appear. For a detailed discussion, see [START_REF] Klein | Numerical Study of oscillatory Regimes in the Kadomtsev-Petviashvili Equation[END_REF]. = 0.1 and t ≤ 0.4. As a reference solution, we consider the solution calculated with the Hochbruck-Ostermann method with N t = 5000 time steps. The dependence of the normalized L 2 -norm of the difference between this reference solution and the numerical solution is shown in Fig. 3.3 in dependence on the time step with a regression analysis and in dependence on the CPU time. Here we can see clearly the phenomenon of order reduction established analytically for parabolic systems by Hochbruck and Ostermann [START_REF] Hochbruck | Exponential Runge-Kutta Methods for semilinear parabolic Problems[END_REF]. In the stiff regime (here up to errors of order 10 -4 ) DCRK does not converge, the Integrating Factor method shows only first order behavior (as predicted in [START_REF] Hochbruck | Exponential Runge-Kutta Methods for semilinear parabolic Problems[END_REF]), and ETD methods perform best. This implies that the Krogstad method is the most economic for the stiff regime of the KP I equation, which gives the precision one is typically interested in in this context. For higher precisions we find a = 1.87 for the Integrating Factor method, a = 4.20 for DCRK, a = 3.75 for Krogstad's ETD scheme, a = 3.90 for the Cox-Matthews scheme, and a = 3.96 for the Hochbruck-Ostermann scheme.

To study empirically the phenomenon of order reduction in exponential integrators, and to observe the transition from a stiff to a non stiff regime we study the ETD schemes in more detail in Fig. 3.4. This is indicated by the fact that ETD schemes are only of order three in this stiff region instead of order four. It appears that all schemes show a slight order reduction though this is not the case for the Hochbruck-Ostermann method in the parabolic case.

Numerical Solution of Cauchy Problems for the KP II Equation

Doubly periodic solution of KP II

The computation for the doubly periodic solution to KP II is carried out with 2 8 ×2 8 points for (x, y) ∈ [-5π, 5π]×[-5π, 5π] and t ≤ 1 with the parameters b = 1, λ = 0.15, bλ 2 +d = -1, µ 1 = µ 2 = 0.25, ν 1 = -ν 2 = 0.25269207053125, ω 1 = ω 2 = -1.5429032317052, and ϕ 1,0 = ϕ 2,0 = 0. The decrease of the numerical error is shown in Fig. 3.5 in dependence of N t and in dependence on CPU time. From a linear regression analysis in a double logarithmic plot we can see that all schemes are fourth order: one finds a = 4.45 for the Integrating Factor method, a = 4.33 for DCRK, a = 3.97 for Krogstad's ETD scheme, a = 4 for the Cox-Matthews scheme, a = 3.95 for the Hochbruck-Ostermann scheme and a = 4.05 for the IRK4 scheme. As for the Zaitsev solution, DCRK performs best followed by the ETD schemes. We thus confirm Driscoll's results in [START_REF] Driscoll | A composite Runge-Kutta Method for the spectral Solution of semilinear PDEs[END_REF] on the efficiency of his method for this example. The absence of order reductions indicates again that the exact solution tests the equation in a non-stiff regime. IRK4 is competitive for larger time steps in this case since only very few iterations (1-3) are needed.

Small dispersion limit for KP II

We consider the same initial data and the same methods as for KP I. In Fig. 3.6 the time evolution of these data can be seen. The solution develops tails this time in negative x-direction. The steepening of the wave fronts happens at essentially the same time, but the gradients are stronger in the vicinity of the tails (see [START_REF] Klein | Numerical Study of oscillatory Regimes in the Kadomtsev-Petviashvili Equation[END_REF]). This is also where the stronger oscillations appear. The computation is carried out with 2 11 × 2 9 points for (x, y) ∈ [-5π, 5π] × [-5π, 5π], = 0.1 and t ≤ 0.4. As a reference solution, we consider the solution calculated with the Hochbruck-Ostermann method with N t = 5000 time steps. The dependence of the normalized L 2 -norm of the difference between this reference solution and the numerical solution on N t and on CPU time is shown in Fig. 3.7. We obtain similar results as in the small dispersion limit of KP I: for typical accuracies one is interested in in this context, DCRK does not converge, the Integrating Factor method shows only a first order behavior, and the ETD methods perform best. In the non-stiff regime we find a = 1.06 for the Integrating Factor method, a = 4.42 for DCRK, a = 4.15 for Krogstad's ETD scheme, a = 4.07 for the Cox-Matthews scheme, and a = 3.99 for the Hochbruck-Ostermann scheme. Once again, we study empirically the phenomenon of order reduction in exponential integrators, and observe a transition from a stiff to a non stiff region (Fig. 3.8), indicated by the fact that ETD schemes are only of order three in this stiff region instead of order four.

Davey-Stewartson II Equation

In this section we perform a similar study as for KP of the efficiency of fourth methods in solving Cauchy problems for the DS II equations. Numerical studies of exact solutions of the DS system were done in [START_REF] White | Numerical simulation of solitons and dromions in the Davey-Stewartson system[END_REF] and [START_REF] Mcconnell | Localised coherent solutions of the DSI and DSII equations-a numerical study[END_REF]. In [START_REF] Besse | Numerical Study of the Davey-Stewartson System[END_REF] blow-up for DS II was studied for the analytically known blow-up solution by Ozawa, see Theorem 1.2.

We consider the equations again on T 2 × R. Due to the ellipticity of the operator in the equation for Φ, it can be inverted in Fourier space in standard manner by imposing periodic boundary conditions on Φ as well. As before this case contains Schwartzian functions that are periodic for numerical purposes. Notice that solutions to the DS equations for Schwartzian initial data stay in this space at least for finite time in contrast to the KP case. Using Fourier transformations, Φ can be eliminated from the first equation by a transformation of the second equation in (3.2) and an inverse transformation. With (3.4) we have

Φ = -2F -1 k 2 x k 2 x + k 2 y F |u| 2 , (3.10)
which leads in (3.2) as for KP to a nonlocal equation with a Fourier multiplier. This implies that the DS equation requires an additional computational cost of 2 two-dimensional FFT per intermediate time step, thus doubling the cost with respect to the standard 2d NLS equation. Notice that from a numerical point of view the same applies to the elliptic-elliptic DS equation that is not integrable. Our experiments indicate that except for the additional FFT mentioned above, the numerical treatment of the 2d and higher dimensional NLS is analogous to the DS II case studied here. The restriction to this case is entirely due to the fact that one can hope for an asymptotic description of the small dispersion limit in the integrable case. Thus we study initial data of the form u 0 (x, y) = a(x, y) exp(ib(x, y)/ ) with a, b ∈ R, i.e., the semi-classical limit well known from the Schrödinger equation. Here we discuss only real initial data for convenience. Notice that Sung's condition (1.47) has been established for the DS II equation with = 1. The coordinate change x = x/ , t = t/ transforms the DS equation (3.2) to this standard form. This implies for the initial data u 0 = exp(-x 2 -ηy 2 ) we study for the small dispersion limit of the focusing DS II system in this section that condition (1.47) takes the form

1 2 η ≤ 1 8 √ 5 -1 2 2 ∼ 0.0477.
This condition is not satisfied for the values of and η we use here. Nonetheless we do not observe any indication of blow-up on the shown timescales. One of the reasons might be that the rescaling with above also rescales the critical time for blow-up by a factor 1/ . In addition it is expected that the dispersionless equations will for generic initial data have a gradient catastrophe at some time t c < ∞, and that the dispersion will regularize the solution for small times t > t c > 0. There are no analytic results in this context, however.

Typically, multi-soliton solutions to the DS II equation will be, as in the KP case, localized in one spatial direction and infinitely extended in another, the lump solution (1.53) is localized in two spatial directions, but with an algebraic fall off towards infinity. Thus these are again not convenient to test codes based on Fourier methods as in the KP case. Since the study of the small dispersion limit below indicates that the time steps have to be chosen sufficiently small for accuracy reasons such that in contrast to KP no order reduction is observed, we will not study any exact solutions here.

The DS II equation is discretized first in space by Fourier spectral method which leads to an equation of the form (2.17) with

L = i αk 2 y -k 2 x (3.11) and N (v, t) = 2iρ F φ + |u| 2 .u (3.12)
with Φ given by (6.6).

We consider here the same methods as for KP, and in addition second and fourth order splitting schemes as given in [START_REF] Yoshida | Construction of higher order symplectic integrators[END_REF], since the DS equation can be split into

i ∂ t u = 2 (-∂ xx u + α∂ yy u), ∂ xx Φ + α∂ yy Φ + 2 |u| 2 xx = 0 (3.13) i ∂ t u = -2ρ Φ + |u| 2 u, (3.14)
which are explicitly integrable, the first two in Fourier space, equation (3.14) in physical space since |u| 2 is a constant in time for this equation.

We mainly include a second order scheme here because of the additional computational cost due to the function Φ in the DS system. This could make a second order scheme competitive in terms of CPU time because of the lower number of FFT used per time step. We will see in the next section ( see Fig. 3.10) that this is not the case.

Small Dispersion Limit for DS II in the defocusing Case

We consider initial data u 0 of the form

u 0 (x, y) = e -R 2 , where R = x 2 + ηy 2 with η = 1 (3.15)
The defocusing effect of the defocusing DS II equation for these initial data can be seen in Fig. 3.9, where |u| 2 is shown for several values of t. The compression of the initial pulse into some almost pyramidal shape leads to a steepening on the 4 sides parallel to the coordinate axes and to oscillations in these regions. The computations are carried out with 2 10 × 2 10 points for (x, y) ∈ [-5π, 5π] × [-5π, 5π], = 0.1 and t ≤ 0.8. To determine a reference solution, we compute solutions with 6000 time steps with the ETD, the DCRK and the IF schemes and take the arithmetic mean. The dependence of the normalized L 2 norm of the difference of the numerical solutions with respect to this reference solution on N t and on CPU time is shown in Fig. 3.10. A linear regression shows that all fourth order schemes show a fourth order behavior except for IRK4 (a = 1.78), as is obvious from the

Small Dispersion Limit for the focusing DS II Equation

For the focusing DS II in the small dispersion limit we consider initial data of the form (3.15) with η = 0.1 and the same methods as before. The focusing effect of the equation can be clearly recognized in Fig. 3.11. The initial peak grows until a breakup into a pattern of smaller peaks occurs. It is crucial to provide sufficient spatial resolution for the central peak. As for the 1+1dimensional focusing NLS discussed in [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger Equation[END_REF], the modulational instability of the focusing DS II leads to numerical problems if there is no sufficient resolution for the maximum. In [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger Equation[END_REF] a resolution of 2 13 modes was necessary for initial data e -x 2 and = 0.1 for the focusing NLS in 1 + 1 dimensions. The possibility of blow-up in DS requires at least the same resolution despite some regularizing effect of the nonlocality Φ. With the computers we could access, a systematic study of time integration schemes with a resolution of 2 13 × 2 13 was not possible in Matlab. Thus we settled for initial data close to the one-dimensional case, which allowed for a lower resolution, see the next section for the Fourier coefficients. The computation is carried out with 2 12 × 2 11 points for (x, y) ∈ [-5π, 5π] × [-5π, 5π], = 0.1 and t ≤ 0.6. To determine a reference solution, we compute solutions with 6000 time steps with the ETD, the DCRK and the IF schemes and take the arithmetic mean. The dependence of the normalized L 2 norm of the difference of the numerical solutions with respect to this reference solution on N t and on CPU time is shown in Fig. 3.12. All schemes except the time splitting scheme show a fourth order behavior, as can be seen from the straight lines with slope a = 3.36 for the Integrating Factor method, a = 3.93 for DCRK, a = 4.06 for Krogstad's ETD scheme, a = 4.11 for the Cox-Matthews scheme, a = 4.13 for the Hochbruck-Ostermann scheme, and a = 2.5 for the fourth order time splitting method. We conclude that in this context DCRK performs best, followed by the ETD schemes. We do not present results for the IRK4 scheme here since it was computationally too expensive.

An Indicator of Accuracy: The Numerical Conservation of the L 2 -Norm

The complete integrability of the KP and the DS equations implies the existence of many or infinitely many conserved quantities (see Chap. 1). It can be easily checked that the L 1 -norm and the L 2 -norm of the solution are conserved as well as the energy. We do not use here symplectic integrators that take advantage of the Hamiltonian structure of the equations. Such integrators of fourth order will be always implicit (as discussed in Chap.

2) which will be in general computationally too expensive for the studied equations as the experiment with the implicit IRK4 scheme showed. Moreover it was shown in [START_REF] Berland | Conservation of phase space properties using exponential integrators on the cubic Schrödinger equation[END_REF] that fourth order exponential integrators clearly outperform second order symplectic integrators for the NLS equation.

The fact that the conservation of L 2 -norm and energy is not implemented in the code allows to use the 'numerical conservation' of these quantities during the computation or the lack thereof to test the quality of the code. We will study in this section for the previous examples to which extent this leads to a quantitative indicator of numerical errors. Note that due to the non-locality of the studied PDEs (1.16) and (3.2), the energies both for KP,

E[u(t)] := 1 2 T 2 ∂ x u(t, x, y)) 2 -λ(∂ -1 x ∂ y u(t, x, y)) 2 -2 2 u 3 (t, x, y) dxdy,
and for DS II,

E[u(t)] := 1 2 T 2 2 |∂ x u(t, x, y)| 2 -2 |∂ y u(t, x, y)| 2 -ρ |u(t, x, y)| 4 - 1 2 Φ(t, x, y) 2 + (∂ -1 x ∂ y Φ(t, x, y)) 2 dxdy,
contain anti-derivatives with respect to x. Since the latter are computed with Fourier methods, i.e., via division by k x in Fourier space, this computation is in itself numerically problematic and could indicate problems not present in the numerical solution of 94 the Cauchy problem. Therefore we trace here only the L 2 -norm T 2 |u(t, x, y)| 2 dxdy, where these problems do not appear. In the plots we show the variable test defined as test = |M (t)/M (0) -1|, where M (t) is the numerically computed L 2 -norm in dependence of time.

-3.2 -3 -2.8 -2.6 -2.4 -2.2 -2 -1.8 -1.6 -8 -7 -6 -5 -4 -3 -2 -1 0 -log 10 (Nt) log10(Δ
Notice that numerical conservation of the L 2 -norm can be only taken as an indication of the quality of the numerics if there is sufficient spatial resolution. Therefore we will always present the Fourier coefficients for the final time step for the considered examples. Notice that no dealiasing techniques and no filtering are used. We will discuss below the results for the small dispersion limit. For the KP I example of Fig. 3.2 we get the Fourier coefficients at the final time and the mass conservation shown in Fig. 3.13. It can be seen that the Fourier coefficients decrease in k x -direction to almost machine precision, whereas this is not fully achieved in k y -direction. This is partly due to the necessity to allow extensive studies of the dependence on the time-stepping in finite computing time and thus to keep the spatial resolution low, and partly due to a Gibbs phenomenon mainly in k y -direction due to the formation of the algebraic tails in Fig. 3.2. Mass conservation can be seen to be a viable indicator of the numerical accuracy by comparing with Fig. 3.3: in the range of accuracy in which one is typically interested (∼ 10 -4 ), mass conservation overestimates the actual accuracy by roughly 2 orders of magnitude. It can be seen that it shows also at least a fourth order decrease. The situation is very similar for the small dispersion example for KP II of Fig. 3.6 as can be seen in Fig. 3.14.

For the defocusing DS II equation and the example shown in Fig. 3.9, the Fourier coefficients decrease to machine precision despite the lower resolution than for KP, see Fig. 3.10. One reason for this is the absence of algebraic tails in the solution. The mass shows as for KP at least fourth order dependence on the time step and overestimates the numerical precision by roughly two orders of magnitude. This is not true for the splitting scheme for which mass conservation is no indication of the numerical precision at all. This seems to be due to the exact integration of the equations (3.14) into which DS is split (for one of them the L 2 norm is constant). The found numerical mass does not appear to reflect the splitting error that is the reason for the numerical error here. For the small dispersion example for the focusing DS II equation of Fig. 3.9 it can be seen in Fig. 3.16 that spatial resolution is almost achieved. There is a certain lack of resolution in the k x direction which leads to the formation of some structure close to k y = 0. This is related to the modulational instability of solutions to the focusing DS II equation. It will disappear for higher resolutions. Numerical conservation of the L 2 -norm of the solution overestimates numerical accuracy by 2-3 orders of magnitude for an error of the order of 10 -3 . Once more it cannot be used as an indicator for the numerical error in the splitting case, where it is almost independent of the time step. For the other cases numerical conservation of the L 2 -norm shows a dependence on N t between fourth and fifth order. This indicates as for the NLS case in [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger Equation[END_REF] that the numerical error has a divergence structure which leads to a higher order decrease of the L 2 -norm than for the actual error. This 96 behavior is also present in the above examples, but less pronounced.

97

Conclusion

It was shown in this chapter that fourth order time stepping schemes can be efficiently used for higher dimensional generalizations of the KdV and the NLS equations, where the 99 stiffness of the system of ODEs obtained after spatial discretization can be a problem. Implicit schemes as IRK4 are computationally too expensive in the stiff regime, whereas standard explicit schemes as RK require for stability reasons too restrictive requirements on the time steps for the KP and DS equations. IMEX schemes do not converge in general for similar reasons, due to the non-localities in the PDEs which lead to singular Fourier multipliers. Driscoll's composite RK variant is generally very efficient if the studied system is not too stiff, but fails to converge for strong stiffness. Exponential integrators do not have this problem. The order reduction phenomenon is a considerable problem for IF schemes in the stiff regime, but less so for ETD schemes. The Hochbruck-Ostermann method performs in general best, but the additional stage it requires is in practice not worth the effort in comparison with Krogstad's or Cox-Matthews' method. The computation of the φ-functions in ETD is inexpensive for the studied problems since it has to be done only once. Since stiffness is not the limiting factor for DS II, all schemes perform well in this context. But the modulational instability of the focusing case requires high spatial resolution we could not achieve in Matlab on the used computers for more general initial data. Thus the code has been parallelized to allow the use of higher spatial resolution without allocating too much memory per processor. This will be discussed in the next chapter.

Chapter 4

Parallel Computing

Although the performance of sequential computers increases incredibly fast, it is insufficient for a large number of challenging applications, in which for example, problems are so large and/or complex that it is impractical or impossible to solve them on a single computer, especially given limited computer memory. Applications requiring much higher performance include numerical simulations in industry and research as well as commercial applications. Parallel programming does not only offer the possibility to solve such problems, but to save also time and money. Indeed, throwing more resources at a task will shorten the time it takes to complete, with potential cost savings, since making a single processor faster is increasingly expensive compared to building parallel clusters from cheap, commodity components. In this short chapter, we introduce some basic notions and terminology related to parallel computing and describe the parallel implementation of our codes.

From Sequential to Parallel Computation

Generally, software is written for sequential computations, in which a code is running on a single computer, having a single Central Processing Unit (CPU). A problem is thus broken into a discrete series of instructions, which are executed one after another, making only one instruction executable at any moment in time. Roughly speaking, parallel computing is the simultaneous use of multiple computational resources to solve a computational problem, which can be a single computer with multiple processors, or different computers connected by a network, or even a combination of both. For an illustration, see Fig. 4.2, in which we want to work in parallel computing on 4 processes, it can be achieved by using for example a single computer with a processor "Quad cores" (case 1), or by using four mono-process computers, linked by a network (case 2) or also by two computers with "Dual cores" processors (case 3). In parallel computations, a problem is broken into discrete parts that can be solved concurrently, each part is then broken down to a series of instructions, and instructions from each part execute simultaneously on different CPUs. Thus, a parallel computer (or multiprocessor system) is a computer utilizing more than one processor.

Parallel computing has its own terminology, and we give here few of it, that will be used in the following. A task is typically a program or program-like set of instructions that is executed by a processor; a parallel task can be executed by multiple processors safely. The scalability of a code is the ability of a parallel code to demonstrate a good scaling. The scaling shows how the time of completion of a problem varies with the number of processors for a fixed total problem size. It is one of the simplest and most widely used indicators to know if a parallel program is performant.

We distinguish two main different architectures: shared and distributed memory architectures.

In a distributed memory architecture, the system is composed of single-processor nodes with local memory. Since processors have their own local memory, memory addresses in Figure 4.3: Parallel Computation one processor do not map to another processor, and so there is no concept of global address space across all processors. Therefore, these systems require a communication network to connect inter-processor memory.

The most important characteristic of this architecture is that access to the local memory is faster than to remote memory. It is nevertheless the challenge for the programmer to assign data to the processors such that most of the data accessed during the computation are already in the node's local memory; and when a processor needs access to data in another processor, it is usually the task of the programmer to explicitly define how and when data is communicated.

The major advantage of distributed memory systems is their ability to scale to a very large number of nodes. In addition, each processor can rapidly access its own memory without interference. Its main disadvantage is that such systems are very difficult to program, the programmer is responsible for many of the details associated with data communication between processors.

In contrast, a shared memory architecture provides (in hardware) a global address space, i.e., all memory locations can be accessed by all processors, via usual load and store operations. Access to a remote location results in a copy of the appropriate cache line in the processor's cache. Consequently, such a system is much easier to program. However, shared memory systems can only be scaled to moderate numbers of processors, because adding more CPUs can geometrically increase traffic on the shared memory-CPU path, and for cache coherent systems, geometrically increase traffic associated with cache/memory management. In addition these systems are expensive. It becomes increasingly difficult and expensive to design and produce shared memory machines with ever increasing num-bers of processors. We work on distributed memory in the following.

Parallelization of the Code

Applications that scale to a large number of processors usually perform computations on large data domains. A prerequisite for parallel numerical algorithms is that sufficient independent computations can be identified for each processor, that require only small amounts of data to be communicated between periods of independent computations.

Identify Sufficient Independent Computations for each Processor: Data Parallel Programming

The basic operations on each object in the data domain (vector, matrix...) have to be executed safely in parallel by the available processors. The scheduling of operations to processors is determined by a domain decomposition specified by the programmer. Processors then execute those operations independently of each others, and determine new values for elements stored in their own local memory. While processors execute an operation, they may need values from other processors. The domain decomposition has thus to be chosen such that the distribution of operations is balanced and the communication is minimized.

The programming model for a distributed memory architecture is the Message Passing Model. Typically in this model, the domain decomposition is implemented by developing a code describing the local computations and local data structures of a single process. Thus, global arrays have to be divided and each parallel task then works on a portion of the data. This handling of global data structures is called data distribution.

There are different ways to partition data, we choose to proceed in the following way: Denoting by

x n = 2πnLx Nx , y m = 2πmLy
Ny , the respective discretizations of x and y in the corresponding computational domain x × y ∈ [-L x π, L x π] × [-L y π, L y π], u is then represented by a N x × N y matrix, which is distributed among processors as follows: Let P 0 , P 1 , ..., P np-1 be the n p processors we can access, then each processor P i will receive N x × Ny np elements of u corresponding to the elements

u(1 : N x , (i -1). N y n p + 1 : i. N y n p ) (4.1)
in the global array. For example, for N x = N y = 8, n p = 4, we distribute the data as shown in Fig. 4.4.

The message passing model is based on a set of processes with private data structures.

Processes communicate by exchanging messages with special send and receive operations. Notice that it can also be used on shared memory computers.

Computations on the global arrays also have to be transformed, e.g., by adapting the loop A variety of message passing libraries have been available since the 1980s, but these implementations differed substantially from each other making it difficult for programmers to develop portable applications.

In 1992, the MPI Forum was formed and became rapidly the standard choice for message passing parallel code implementations.

Minimize Communications/Exchange of Data between Processors

As explained before, in distributed parallel programming, different processors work on completely independent data. Consequently, processors have to communicate with each other to exchange data by Message Passing Interface (MPI). Communications are done explicitly by calling sub-routines of the MPI library.

The only part of our code that requires communications is the computation of the discrete two-dimensional Fourier transforms and its inverse. To perform them, we use the transposition approach. The latter allows to use highly optimized single processor one-dimensional FFT routines, that are normally found in most architectures, and a transposition algorithm can be easily ported to different distributed memory architectures. We use the well known FFTW library because its implementation is close to optimal for serial FFT computation, see [43]. Roughly speaking, a 2d FFT does 1-d FFTs on all rows and then on all columns of the initial global array.

We thus first transform in x direction, each processor transforms all the dimensions of the data that are completely local to it by using the 1-d FFTW routine; and then, because our subroutine FFTW does only one direction, in order to evaluate the Fourier Transform in the y direction, one needs to transpose the global matrix, thus interchanging the directions.

To do this, we use MPI ALLTOALL communications to perform the data transposition.

Essentially, the way we do this is to do an in-place transpose of the data local to each processor (itself a non-trivial problem), exchange blocks of data with every other processor, and then complete the transpose by a further (precomputed) permutation on the data. For example, consider that we have N x = N y = n p = 4, we then obtain the transpose of the initial matrix distributed correctly amongst processors. After the transpose, the Fourier transform û is distributed on the processes so that process P i holds the elements corresponding to û((i -1)

N x n p + 1 : iN x n p , 1 : N y),
on which we are now able to call again the 1d FFTW subroutine, and we have computed the 2d-FFT (one could also transpose back, but this is optional; it is not done here to avoid additional communications).

A bottleneck experienced in scaling pseudospectral codes to large processor counts is in performing output to disk. This often occurs significantly slower than computation. To alleviate this, we use the simplest method that we know, that is to produce independently one binary file per process per saved data. One of possible disadvantage of this method is the necessity to post-treat the data after the computation to visualize the results, this is done by a script already implemented, and requires only few minutes of computation more.

We show below in Fig. 4.6 the scaling of our code which is better as the number of modes increases.

Numerical Study of Blowup in the Davey-Stewartson System

There exist many explicit solutions for the integrable cases of the DS system which thus allow to address the question about the long time behavior of solutions for given initial data.

For the famous Korteweg-de Vries (KdV) equation, it is known that general initial data are either radiated away or asymptotically decompose into solitons. The DS II system and the two-dimensional integrable generalization of KdV known as the Kadomtsev-Petviashvili I (KP I) equation have so-called lump solutions, a two-dimensional soliton which is localized in all spatial directions with an algebraic fall off towards infinity. For KP I it was shown [START_REF] Ablowitz | On the inverse scattering and direct linearizing transforms for the Kadomtsev-Petviashvili equation[END_REF] that small initial data asymptotically decompose into radiation and lumps. It is conjectured that this is also true for general initial data.

For the defocusing DS II global existence in time was shown by Fokas and Sung [START_REF] Fokas | The inverse spectral method for the KP I equation without the zero mass constraint[END_REF] for solutions of certain classes of Cauchy problems. These initial data will simply disperse. The situation is more involved for the focusing case. Pelinovski and Sulem [START_REF] Pelinovsky | Eigenfunctions and Eigenvalues for a scalar Riemann-Hilbert Problem associated to inverse Scattering[END_REF] showed that the lump solution is spectrally unstable. In addition the focusing NLS equations in 2 + 1 dimensions have the critical dimension, i.e., solutions from smooth initial data can have blowup. This means that the solutions lose after finite time the regularity of the initial data, a norm of the solution or one of its derivatives becomes infinite. For focusing NLS equations in 2 + 1 dimensions, it is known that blowup is possible if the energy of the initial data is greater than the energy of the ground state solution, see e.g. [START_REF] Sulem | The nonlinear Schrödinger equation[END_REF] and references therein, and [START_REF] Merle | On universality of blow-up profile for l 2 critical nonlinear Schrödinger equation[END_REF] for an asymptotic description of the blowup profile. For the focusing DS II equation Sung [START_REF] Sung | Long-time decay of the solutions of the Davey-Stewartson II equations[END_REF] gave a smallness condition on the Fourier transform F[u 0 ] of the initial data to establish global existence in time for solutions to Cauchy problems given by the inequality 1.47. It is not known whether there is generic blowup for initial data not satisfying this condition, nor whether the condition is optimal. Since the initial data studied in this chapter are not in this class, we cannot provide further insight into this question. An explicit solution with blowup for lump-like initial data was given by Ozawa [START_REF] Ozawa | Exact Blow-up Solutions to the Cauchy Problem for the Davey-Stewartson Systems[END_REF]. It has an L ∞ blowup in one point (x c , y c , t c ) and is analytic for all other values of (x, y, t). It is unknown whether this is the typical blowup behavior for the focusing DS II equation.

From the point of view of applications, a blowup of a solution does not mean that the studied equation is not relevant in this context, but just indicates the limit of the used approximation. It is thus of particular interest, not only in mathematics, but also in physics, since it shows the limits of the applicability of the studied model. This breakdown of the model will also in general indicate how to amend the used approximations.

In view of the open analytical questions concerning blowup in DS II solutions, we study the issue in the present chapter numerically, which is a highly non-trivial problem for several reasons: first DS is a purely dispersive equation which means that the introduction of numerical dissipation has to be avoided as much as possible to preserve dispersive effects such as rapid oscillations. This makes the use of spectral methods attractive since they are known for minimal numerical dissipation and for their excellent approximation properties for smooth functions. But the algebraic falloff of both the lump and the Ozawa solution leads to strong Gibbs phenomena at the boundaries of the computational domain if the solutions are periodically continued there. We will nonetheless use Fourier spectral methods because they also allow for efficient time integration algorithms which should be ideally of high order to avoid a pollution of the Fourier coefficients due to numerical errors in the time integration. An additional problem is the modulational instability of the focusing DS II equation, i.e., a self-induced amplitude modulation of a continuous wave propagating in a nonlinear medium, with subsequent generation of localized structures, see for instance [START_REF]Nonlinear fiber optics[END_REF][START_REF] Cross | Pattern formation outside of equilibrium[END_REF][START_REF] Forest | Geometry and modulation theory for the periodic nonlinear Schrödinger equation[END_REF] for the NLS equation. Thus to address numerically questions of stability and blowup of its solutions, high resolution is needed which cannot be achieved on single processor computers. Therefore we use parallel computers1 to study the related questions. The use of Fourier spectral method is also very convenient in this context, since for a parallel spectral code only existing optimized serial FFT algorithms are necessary.

In addition such codes are not memory intensive, in contrast to other approaches such as finite difference or finite element methods. The first numerical studies of DS were done by White and Weideman [START_REF] White | Numerical simulation of solitons and dromions in the Davey-Stewartson system[END_REF] using Fourier spectral methods for the spatial coordinates and a second order time splitting scheme. Besse, Mauser and Stimming [START_REF] Besse | Numerical Study of the Davey-Stewartson System[END_REF] used essentially a parallel version of this code to study the Ozawa solution and blowup in the focusing elliptic-elliptic DS equation. McConnell, Fokas and Pelloni [START_REF] Mcconnell | Localised coherent solutions of the DSI and DSII equations-a numerical study[END_REF] used Weideman's code to study numerically DS I and DS II, but did not have enough resolution to get conclusive results for the blowup in perturbations of the lump in the focusing DS II case. In this chapter we repeat some of their computations with considerably higher resolution.

We use a parallelized version of a fourth order time splitting scheme which was studied for DS in the third chapter, this method being very convenient for parallel computing, because of easy coding (loops) and low memory requirements. It is applicable to the DS II equation, which can be split into

i∂ t u = (-∂ xx u + α∂ yy u), ∂ xx Φ + α∂ yy Φ + 2∂ xx |u| 2 = 0, (5.1
)

i∂ t u = -2ρ Φ + |u| 2 u, ( 5.2)
which are explicitly integrable, the first two in Fourier space, equation (5.2) in physical space since |u| 2 is a constant in time for this equation.

Notice that the splitting method in the form (5.2) conserves the L 2 norm: the first equation implies that its solution in Fourier space is just the initial condition (from the last time step) multiplied by a factor e iφ with φ ∈ R. Thus the L 2 norm is constant for solutions to this equation because of Parseval' s theorem. The second equation as mentioned conserves the L 2 norm exactly. Thus the used splitting scheme has the conservation of the L 2 norm implemented. As we will show in the following, this does not guarantee the accuracy of the numerical solution since other conserved quantities as the energy the conservation of which is not implemented might not be numerically conserved. In fact we will show that the numerically computed energy provides a valid indicator of the quality of the numerics.

Obviously it is non-trivial to decide numerically whether a solution blows up or whether it just has a strong maximum. To allow to make nonetheless reliable statements, we perform a series of tests for the numerics. First we test the code on known exact solutions with algebraic falloff, the lump and the Ozawa solution. We establish that energy conservation can be used to judge the quality of the numerics. It is shown that the splitting code continues to run in general beyond a potential blowup which makes difficult to decide whether there is blowup or just a strong maximum of the solution. We argue at examples for the quintic NLS in 1 + 1 dimensions (which is known to have blowup solutions) and the Ozawa solution that energy conservation is a reliable indicator in this case since the energy of the solution changes completely after a blowup, whereas it will be in accordance with the numerical accuracy after a strong maximum. Thus we reproduce well known blowup cases in this way and establish with the energy conservation a criterion to ensure the accuracy of the numerics also in unknown cases. Then we study perturbations of the lump and the Ozawa solution to see when blowup is actually observed.

Lump solution of the focusing DS II equation

To test the performance of the code, we first propagate initial data from known exact solutions and compare the numerical and the exact solution at a later time.

The focusing DS II equation has solitonic solutions which are regular for all x, y, t, and which are localized with an algebraic falloff towards infinity, known as lumps [START_REF] Arkadiev | Inverse Scattering Transform Method and Soliton Solutions for Davey-Stewartson II Equation[END_REF]. The single lump is given by

u(x, y, t) = 2c exp -2i(ξx -ηy + 2(ξ 2 -η 2 )t) |x + 4ξt + i(y + 4ηt) + z 0 | 2 + |c| 2 (5.3)
where (c, z 0 ) ∈ C 2 and (ξ, η) ∈ R 2 are constants. The lump moves with constant velocity (-2ξ, -2η) and decays as (x 2 + y 2 ) -1 for x, y → ∞.

Blowup for the quintic NLS in 1 + 1 dimensions and the focusing DS II

It is known that focusing NLS equations can have solutions with blowup, if the nonlinearity exceeds a critical value depending on the spatial dimension. For the 1+1 dimensional case, the critical nonlinearity is quintic, for the 2 + 1 dimensional it is cubic, see for instance [START_REF] Sulem | The nonlinear Schrödinger equation[END_REF] and references therein. Thus the focusing DS II equation can have solutions with blowup.

In this section we will first study numerically blowup for the 1 + 1 dimensional quintic NLS equation, and then numerically evolve initial data for a known exact blowup solution to the focusing DS II equation due to Ozawa [START_REF] Ozawa | Exact Blow-up Solutions to the Cauchy Problem for the Davey-Stewartson Systems[END_REF]. We discuss some peculiarities of the fourth order splitting scheme in this context.

Blowup for the quintic one-dimensional NLS

The focusing quintic NLS in 1 + 1 dimensions has the form

i∂ t u + ∂ xx u + |u| 4 u = 0, (5.4)
where u ∈ C depends on x and t (we consider again solutions periodic in x). This equation is not completely integrable, but assuming the solution is in L 2 , has conserved L 2 norm and, provided the solution u ∈ H 2 , a conserved energy,

E[u] = R 1 2 |∂ x u| 2 - 1 6 |u| 6 dx. (5.5)
It is known that initial data with negative energy blow up for this equation in finite time, and that the behavior close to blowup is given in terms of a solution to an ODE, see [START_REF] Merle | On universality of blow-up profile for l 2 critical nonlinear Schrödinger equation[END_REF].

As discussed in sect. 2.1, the splitting scheme we are using here has the property that the L 2 norm is conserved. Thus the quality of the numerical conservation of the L 2 norm gives no indication on the accuracy of the numerical solution. However as discussed in [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger Equation[END_REF], conservation of the numerically computed energy gives a valid criterion for the quality of the numerics: in general it overestimates the L ∞ numerical error by two orders of magnitude at the typically aimed at precisions. If we consider as in [START_REF] Stinis | Numerical computation of solutions of the critical nonlinear Schrödinger equation after the singularity[END_REF] for the quintic NLS the initial data u 0 (x) = 1.8i exp(-x 2 ), the energy is negative. We compute the solution with L x = 5 and N x = 2 15 with N t = 10 4 time steps. The result can be seen in Fig. 5.3 (to obtain more structure in the solution after the blow up due to a less pronounced maximum, the plot on the left was generated with the lower spatial resolution N = 2 12 ). The initial data clearly get focused to a strong maximum, but the code does not break. We note that this is in contrast to other fourth order schemes tested for 1 + 1 dimensional NLS equations in [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger Equation[END_REF], which typically produce an overflow close to the blowup. But clearly the solution shows spurious oscillations after the time t c ∼ 0.155. In fact the numerically computed energy, which will always be time-dependent due to unavoidable numerical errors, will be completely changed after this time. We consider

∆ E = 1 - E(t) E(0) , ( 5.6)
where E(t) is the numerically computed energy (5.5) and get for the example in Fig. 5.3 the behavior shown in Fig. 5.4. At the presumed blowup at t c ∼ 0.155 as in [START_REF] Stinis | Numerical computation of solutions of the critical nonlinear Schrödinger equation after the singularity[END_REF], the energy jumps to a completely different value. Thus this jump can and will be used to indicate blowup. To illustrate the effects of a lower resolution in time and space imposed by hardware limitations for the DS computations, we show this quantity for several resolutions in Fig. 5.4. If a lower resolution in time is used as in some of the DS examples in this paper, the jump is slightly smoothed out. But the plateau is still reached at essentially the same time which indicates blowup. Thus a lack of resolution in time in the given limits will not be an obstacle to identify a possible singularity. The reason for this is the use of a fourth order scheme that allows to take larger time steps. We will present computations with different resolutions to illustrate the steepening of the energy jump as above if this is within the limitations imposed by the hardware. We show the modulus of the Fourier coefficients for N = 2 12 and N = 2 15 before and after the critical time in Fig. 5.5. It can be seen that the solution is well resolved before blowup in the latter case, and that the singularity leads to oscillations in the Fourier coefficients. A lack of spatial resolution as for N = 2 12 in Fig. 5.5 triggers the modulation instability close to the blowup and at later times as can be seen from the modulus of the Fourier coefficients that increase for larger wavenumbers. Therefore we always aim at a sufficient resolution in space even for times close to a blowup. After this time the modulation instability will be present in the spurious solution produced by the splitting scheme as we will show for an example. Remark 5. 1 Stinis [89] has recently computed singular solutions to the focusing quintic nonlinear Schrödinger equation in 1 + 1 dimensions. This equation has solutions in L ∞ L 2 that may not be unique for given smooth initial data and that may exhibit blowup of the L ∞ H 1 norm. Following Tao [START_REF] Tao | Global existence and uniqueness results for weak solutions of the focusing mass-critical nonlinear Schrödinger equation[END_REF], Stinis [START_REF] Stinis | Numerical computation of solutions of the critical nonlinear Schrödinger equation after the singularity[END_REF] has used a selection criteria to pick a solution after the blow up time of the L ∞ H 1 norm. They suggest that 'mass' is ejected (which means that the L 2 norm is changed) at times where the L ∞ H 1 norm blows up. The splitting scheme studied here in contrast produces a weak solution with a different energy since the L 2 norm conservation is built in.

Blowup in the Ozawa solution

For the focusing DS II equation, an exact solution was given by Ozawa [START_REF] Ozawa | Exact Blow-up Solutions to the Cauchy Problem for the Davey-Stewartson Systems[END_REF] which is in L 2 for all times with an L ∞ blowup in finite time. His results are recalled in the Theorem 1.2.

We thus consider initial data of the form

u(x, y, 0) = 2 exp -i(x 2 -y 2 ) 1 + x 2 + y 2 (5.

Perturbations of the lump solution

In this section we consider perturbations of the lump solution (5.3). First we propagate initial data obtained from the lump after multiplication with some scalar factor. Then we consider a perturbation with a Gaussian and a deformed lump.

Perturbation of the lump by rescaled initial data

We first consider rescaled initial data from the lump (5. They again decrease by almost 6 orders of magnitude. To illustrate the modulational instability at a concrete example, we show the Fourier coefficients after the critical time in Fig. 5.13. It can be seen that the modulus of the coefficients of the high wavenumbers increases instead of decreasing as to be expected for smooth functions. This indicates once more that the computed solution after the blowup time has to be taken with a grain of salt.

For A = 0.9, the initial pulse travels in the same direction as the exact solution, but loses speed and height and is broadened, see Fig. 5.14. It appears that this modified lump just disperses asymptotically. The solution can be seen in Fig. 5.15. Its Fourier coefficients in Fig. 5.16 show that the resolution of the initial data is almost maintained.

Perturbation of the lump with a Gaussian

We consider an initial condition of the form u(x, y, -6) = u l + B exp(-(x 2 + y 2 )), B ∈ R.

(5.8)

For B = 0.1 and N t = 1000, we show the solution at different times in Fig. 5.17. The solution travels at the same speed as before, but its amplitude varies, growing and decreasing successively, see Fig. 5.18. The time evolution of the energy can be seen in Fig. 5.18.

There is no indication of blowup in this example. The solution appears to disperse for t → ∞. The Fourier coefficients at t = 6 in Fig. 5.19 show the wanted spatial resolution.

A similar behavior is observed if a larger value for the amplitude of the perturbation is chosen, e.g., B = 0.5, see Fig.

Deformation of the Lump

We consider initial data of the form u(x, y, -6) = u l (x, κy, -6), (5.9) For κ = 0.9, the resulting solution loses speed and width as can be seen in Fig. 5.21. Its height and energy grow, but both stay finite, see Fig. 5.22. It is possible that the solution eventually blows up, but not on the time scales studied here. The Fourier coefficients at t = 0 in Fig. 5.23 show the wanted spatial resolution.

For κ = 1.1, we observe the opposite behavior in Fig. 5.24. The solution travels with higher speed than the initial lump and is broadened. The energy does not show any sudden change, see Fig.

Perturbations of the Ozawa solution

In this section we study as for the lump in the previous section various perturbations of initial data for the Ozawa solution to test whether blowup is generic for the focusing DS II equation.

Perturbation of the Ozawa solution by multiplication with a scalar factor

We consider initial data of the form For C = 0.9, the initial pulse grows until it reaches its maximal height at t = 0.2501, but there is no indication for blowup, see Fig. Thus for initial data given by the Ozawa solution multiplied with a factor C, we find that for C > 1, blow up seems to occur before the critical time of the Ozawa solution, and for C < 1 the solution grows until t = 0.25 but does not blow up. Consequently the Ozawa initial data seem to be critical in this sense that data of this form with smaller norm do not blow up.

u(x, y, 0) = 2C exp -i(x 2 -y 2 ) 1 + x 2 + y 2 , ( 5

Perturbation of the Ozawa solution with a Gaussian

We consider an initial condition of the form in the Ozawa solution is clearly a generic feature at least for initial data close to Ozawa for the focusing DS II equation.

Deformation of the Ozawa solution

We study deformations of Ozawa initial data of the form For ν = 0.9, we observe a maximum of the solution at t = 0.2441, see Fig. 5.37, followed

Conclusion

In this chapter we have numerically studied long time behavior and stability of exact solutions to the focusing DS II equation with an algebraic falloff towards infinity. We have shown that the necessary resolution can be achieved with a parallelized version of a spectral code. The spatial resolution as seen at the Fourier coefficients was always well beyond typical plotting accuracies of the order of 10 -3 . For the time integration we used an unconditionally stable fourth order splitting scheme. As argued in [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger Equation[END_REF][START_REF] Klein | Fourth order time-stepping for Kadomtsev-Petviashvili and Davey-Stewartson equations[END_REF], the numerically computed energy of the solution gives a valid indicator of the accuracy for sufficient spatial resolution. To ensure the latter, we always presented the Fourier coefficients of the solution at a time before a singularity appeared. In addition we show here that the numerically computed energy indicates blowup by jumping to a different value in cases where the code runs beyond a singularity in time.

After testing the code for exact solutions, the lump and the blowup solution by Ozawa, we showed that both solutions are critical in the following sense: adding energy to it leads to a blowup for the lump, and an earlier blowup time for the Ozawa solution. For initial data with less energy, no blowup was observed in both cases, the initial data asymptotically just seem to be dispersed. This is in accordance with the conjecture in [START_REF] Mcconnell | Localised coherent solutions of the DSI and DSII equations-a numerical study[END_REF] that solutions to the focusing DS II equations either blow up or disperse. In particular the lump is unstable against both blowup and dispersion, in contrast to the lump of the KP I equation that appears to be stable, see for instance [START_REF] Pelinovsky | Eigenfunctions and Eigenvalues for a scalar Riemann-Hilbert Problem associated to inverse Scattering[END_REF]. Note that the perturbations we considered here test the nonlinear regime of the PDE for which so far no analytical results appear to be established.

Chapter 6

Study of the focusing Davey-Stewartson II equation in the small dispersion limit

In this chapter, we first present few results for several computations done in the previous chapter, but this time using the DCRK method, which has been identified in the third chapter to be the most efficient time stepping scheme for the study of DS II. The DCRK code has been parallelized in the way explained in chapter IV; the main difference with the time splitting code being that DCRK method requires the use of two arrays of dimensions N x ×N y /n p in more than in the former. This additional memory required has nevertheless not big impact on the time of computation.We then perform a study of the semiclassical limit of the focusing DS II equation.

DCRK versus Splitting Method

It has been shown in the previous chapter that even if splitting methods are very convenient for the implementation of a parallel code, they have the disadvantage that the code continues to run in general beyond a potential blow-up and present some unwanted peculiarities. This is not the case with the DCRK method, which appears here to be more efficient. We observe a breaking of the code at the blow up time, and the computed energy does not suffer as much from numerical errors than with splitting method. In addition Fourier coefficients reach better precision.

For the Ozawa solution, the computation is carried out with N x = N y = 2 15 modes, L x = L y = 20 and N t = 1000. The code breaks at t = 0.2484, see Fig. 6.1, and in this case, Fourier coefficients decrease to machine precision, which implies that Gibbs phenomena do not play a role, see Fig. 6.2. The modulational instability shows up in this plot in the form of an increasing of the high wave numbers. Since their modulus is still below 136 10 -10 this instability does not affect the results here, but it shows that high resolution is crucial to obtain reliable results in this case. We can notice that the precision is however better than one reached with the time splitting scheme. For the lump solution rescaled by a scalar factor, the DCRK method gives the results below.

For A = 1.1, we recover the blow-up time at t ∼ 1.96, and a jump in the computed energy, see Fig. one, as well in term of recovering of blow up time, smoothness in the energy, as in the Fourier coefficients precision, we choose in the following to investigate the semiclassical limit of the focusing DS II equation with the DCRK method.

Numerical Study of the Semiclassical Limit

In this section we perform a study of the focusing DS II equation in the semiclassical limit, that we did not able to do in the third chapter on a single computer.

We consider the focusing DS II equation in the semiclassical limit ( → 0), and consider initial data of the form u 0 (x, y) = e -R 2 , with R = x 2 + ηy 2 (6.

2)

The computation are carried out with 2 15 × 2 15 points for

x × y ∈ [-5π, 5π] × [-5π, 5π]
for different values of η, and = 0.1. The maximal time of computation is chosen to be t max = 0.6, as twice the apparent break up time of the corresponding dispersionless DS system. The latter is determined by looking for solutions of the form Ψ(x, y, 0, ) = A(x, y)e iS(x,y)/ , with A(x, y) = u(x, y), (6.3) for the DS II equations

i Ψ t + 2 Ψ xx -α 2 Ψ yy + 2ρ Φ + |Ψ| 2 Ψ = 0, Φ xx + βΦ yy + 2 |Ψ| 2 xx = 0, (6.4)
where α = β = 1 and ρ take the values ±1, where 1 is a small dispersion parameter, and where Φ is a mean field.

The dispersionless DS II system ( = 0) is thus given by S t + S 2

x -αS 2 y -2ρ(Φ + u) = 0 u t + 2(S x u x + S xx u) -2α(S y u y + S yy u) = 0 (6.5)

with Φ = -2F -1 k 2 x k 2 x + k 2 y F |u| 2 , ( 6.6)
where F is the 2-dimensional Fourier transform of u.

For η = 0.1, the initial peak grows until its maximal height at t ∼ 0.3. After that smaller humps appear, see Fig. 6.8; as in the case of the one dimensional cubic NLS equation in semiclassical limit. A similar behavior is observed for η = 0.3, see Fig. 6.11.

The behavior of the Fourier coefficients at t max is shown in Fig. 6.10 for the case η = 0.1 and in Fig. 6.13 for η = 0.3. They decrease to machine precision, without showing any increase as before at high wave numbers. For η = 0.5 and η = 0.8, we observe three peaks of the solution, again followed by the apparition of smaller humps, see Fig. 6.14, Fig. 6.15 and 6.16. The behavior of the Fourier coefficients at t max is shown in Fig. 6.17 for the case η = 0.5 and for η = 0.8. They again decrease to machine precision. We compare in Fig. 6.18 the contours of the solutions at the final time of computation for different values of η.

For η = 1, we observe a blow-up at t ∼ 0.28, see Fig. It was shown in this chapter that as already discuss in Chap. III, DCRK method is very well suited for the study of the DS II equation. We showed that parallel computing is an essential tool nowadays for numerical simulations in cases where theories do not yet succeed. The writing of parallelized codes in Fortran 90 allowed to perform a valid numerical study of the focusing DS II equation in the semiclassical limit, leading to the conclusion that a blow up occurs for initial data with radial symmetry, as in the critical The KP equation shows in the small dispersion limit a similar behavior as the Kortewegde Vries equation in the one-dimensional case. We observe dispersive shocks where the corresponding dispersionless equation has shocks. For the defocusing Davey-Stewartson II equation, we again observe a similar behavior as in the one-dimensional case, the solution is broadened, and small oscillations appear before the point of gradient catastrophe.

For the focusing DS II equation, the modulational instability made it necessary to use parallel computing, which enabled the study of the semiclassical limit and also of blow-up phenomena in this case. We worked in collaboration with different cluster centers (as IDRIS, CCRT), and obtained various results. Using the most efficient scheme identified before for this regime (namely Driscoll's composite RK method), we found in the semiclassical limit that dispersive shocks only seem to appear in solutions from initial data without a radial symmetry, whereas blow-up occurs in this case.

On the other hand, we studied the stability with respect to perturbations of the lump solution and the Ozawa solution (which blows up) to the focusing DS II equation and the appearence of blow-up there. This time, we used a time splitting scheme, and some peculiarities of it were observed, as for example, the fact that it produces some artificial continuation of the solution beyond the blow up time. We established that the numerical conservation of a conserved quantity of DS can be used to identify the blow-up time, and that both studied exact solutions are critical in a sense specified in chapter 5. In contrast to the KP I case studied in [START_REF] Ablowitz | On the inverse scattering and direct linearizing transforms for the Kadomtsev-Petviashvili equation[END_REF], where it was shown that small initial data asymptotically decompose into radiation and lumps, we found here that the one-lump of DS II is unstable.

This work gives some indications on important properties of the studied equations, and offers also many perspectives for future work.

Firstly, our work can be improved by establishing scaling laws for various quantities in the vicinity of dispersive shocks. We also studied here the stability of the one-lump soliton for the focusing DS II equation. The same will be done for exact solutions to the NLS equation as solitons and breathers. The blow-up phenomena will be investigated in more detail, too.

Secondly, we studied here only completely integrable systems, whereas few results are known about non-integrable ones. For example the DS system could be studied in the non-integrable cases in more detail numerically with our codes, since they do not use the complete integrability property of them. Here we typically studied the DS system in hyperbolic-elliptic form, which can be written as

i∂ t u + ∂ xx u -∂ yy u = α |u| 2 u + βuΦ (6.7) ∆Φ = ∂ xx |u| 2 (6.8)
for the complex-valued function u(x, y, t) and for the real potential Φ, with α = 1 and β = -2 (which gives the only known integrable hyperbolic-elliptic case, DS II), but values of α and β can be chosen arbitrarily to study non-integrable cases.

analyticity strip is equal to the logarithmic decrement of the Fourier transform at large wavenumbers, that is ln |û(k)| ∼ δk, k → +∞ (6.12) Following these observations, û(k, t) is assumed to be of the form û(k, t) = C(t)k -α(t)k e -δ(t)k . (6.13) This only holds for a function of one spatial variable. However slightly weaker estimates are valid for E(K, t), again for an analytic function u in S α . Thus we assume for a function in several space dimensions that E(K, t) = C(t)K -α(t) e -δ(t)K . (6.14)

One can then trace the temporal behavior of δ(t) in order to obtain evidence for or against blow-up (the problem of blow-up reduces to check if δ(t) vanishes in finite time, which indicates a loss of regularity).

In order to extract δ(t) from direct numerical simulations, a least-square fit is performed on the logarithm of the computed quantity (Fourier transform or energy spectrum), using the functional form: The performance of this method for both dispersionless NLS and DS systems is yet under investigation, and will be the subject of future research, as the other perspectives described above.

1. 1 ω

 1 in dependence of ξ for the 1-d wave equation, with α = 1 . . . . . . 1.2 ω in dependence of ξ for the Airy equation, with α = β = 1 and the convention ω > 0 (blue) on the left, and for the Schrödinger equation with = 1 on the right. . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Evolution of the solution to the linearized gKdV equation, with u 0 = -sech 2 (x). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Time evolution of the solution of (1.25) with u(0, x) = -sech 2 x and n = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Time evolution of the solution of (1.29) with u(0, x) = -sech 2 (x) and ν = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Dispersive shock waves for KdV in small dispersion limit. . . . . . . . 1.7 KdV and asymptotic solution with u 0 = -sech 2 (x) and = 0.01, see [47]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Solution to the defocusing cubic NLS equation in semiclassical limit with u(0, x) = e -x 2 , = 0.5 and t max = 1. . . . . . . . . . . . . . . . . 1.9 Solution to the focusing cubic NLS equation in semiclassical limit with u(0, x) = e -x 2 , = 0.1 and t max = 0.8. . . . . . . . . . . . . . . . 1.10 Time evolution of the solution of (1.32) with u(0, x) = sech (x) and = 0.1, blow-up occurs. . . . . . . . . . . . . . . . . . . . . . . . . . 1.11 Zaitsev solution to the KP I equation, α = 1 and β = 0.5, at t = 0. . 1.12 Doubly periodic hyperelliptic KP II solution of genus 2 with parameters: b = 1, λ = 0.15, bλ 2+ d = -1, µ 1 = µ 2 = 0.25, ν 1 = -ν 2 = 0.25269207053125, ω 1 = ω 2 = -1.5429032317052,and ϕ 1,0 = ϕ 2,0 = 0, at t = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Example of Aliasing effect (x for f 1 , o for f 2 ) on the discrete grid. . . . . 2.2 The precision of the approximation depends on the smoothness of the function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Normalized L 2 -norm of the numerical error in the solution to the KP I equation with initial data given by the Zaitsev solution for several numerical methods, as a function of N t (top) and as a function of CPU time (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Solution to the KP I equation for the initial data u 0 = -∂ x sech 2 (R) where R = √ x 2 + y 2 for several values of t. . . . . . . . . . . . . . . . 3.3 Normalized L 2 -norm of the numerical error for the solution shown in Fig. 3.2 for several numerical methods, as a function of N t (top) and as a function of CPU time (bottom). . . . . . . . . . . . . . . . . . . 3.4 Phenomenon of order reduction in Exponential Integrators for KP I in the small dispersion limit, ETD schemes of Fig. 3.3 . . . . . . . . . 3.5 Normalized L 2 -norm of the numerical error for the time evolution of the doubly periodic solution to the KP II equation for several numerical methods, as a function of N t (top) and as a function of CPU time (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Solution to the KP II equation for the initial data u 0 = -∂ x sech 2 (R) where R = √ x 2 + y 2 for several values of t. . . . . . . . . . . . . . . . 3.7 Normalized L 2 -norm of the numerical error for the solution in Fig. 3.6 for several numerical methods as a function of N t (top) and as a function of CPU time (bottom). . . . . . . . . . . . . . . . . . . . . . 3.8 Phenomenon of order reduction in Exponential Integrators for KPII in the small dispersion limit (see Fig.3.7). . . . . . . . . . . . . . . . 3.9 Solution to the defocusing DS II equation for the initial data u 0 = exp(-R 2 ) where R = √ x 2 + y 2 and = 0.1 for several values of t. . . 3.10 Normalized L 2 -norm of the numerical error for several numerical methods for the situation shown in Fig. 3.9 as a function of N t (top) and of CPU time (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Solution to the focusing DS II equation for the initial data u 0 = exp(-R 2 ) where R = √ x 2 + ηy 2 , η = 0.1 and = 0.1 for several values of t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Normalized L 2 -norm of the numerical error for the example of Fig. 3.11 for several numerical methods as a function of N t (top) and of CPU time (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Non-conservation of the numerically computed L 2 -norm of the solution to the problem considered in Fig. 3.3 in dependence on the time step (top) and the Fourier coefficients for the final time (bottom). . . 3.14 Non-conservation of the numerically computed L 2 -norm of the solution to the problem considered in Fig. 3.7 in dependence on the time step (top) and the Fourier coefficients for the final time (bottom). . . 3.15 Non-conservation of the numerically computed L 2 -norm of the solution to the problem considered in Fig. 3.10 in dependence on the time step (top) and the Fourier coefficients for the final time (bottom). . . 3.16 Non-conservation of the numerically computed L 2 -norm of the solution to the problem considered in Fig. 3.12 in dependence on the time step (top) and the Fourier coefficients for the final time (bottom). . . 4.1 Serial Computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Example of three different parallel architectures to work on 4 different processes, P 0 , P 1 , P 2 and P 3 ; the double arrow representing a network. 4.3 Parallel Computation . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Example of data distribution for N x = N y = 8, n p = 4 . . . . . . . . . 4.5 Reconstructed transposed matrix . . . . . . . . . . . . . . . . . . . . 4.6 Scaling of the code . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Contours of |u| 2 on the left and a plot of ||u exact -u num || 2 on the right in dependence of time for the solution to the focusing DS II equation (3.2) for lump initial data (5.3). . . . . . . . . . . . . . . . . . . . . . 5.2 Fourier coefficients for the situation in Fig. 5.1 at t = 6. . . . . . . . . 5.3 Solution to the focusing quintic NLS (5.4) for the initial data u 0 = 1.8i exp(-x 2 ) with N = 2 12 on the left and N = 2 15 on the right for t > t c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Numerically computed energy for the situation studied in Fig. 5.3 for N = 2 12 on the left and N = 2 15 on the right for several values of N t . At the blowup, the energy jumps. . . . . . . . . . . . . . . . . . . . 5.5 Fourier coefficients for the solution in Fig. 5.3 close to the critical and at a later time for N = 2 12 on the left and N = 2 15 on the right for N t = 10 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Solution to the focusing DS II equation (3.2) for t = 0.075 and t = 0.15 in the first row and t = 0.225 and t = 0.3 below for an initial condition of the form (5.7). . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Time evolution of max(|u num | 2 ) and of u num -u exact 2 for the situation in Fig. 5.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Numerically computed energy E(t) and ∆E = |1 -E(t)/E(0)| (5.6) for the situation in Fig. 5.6. . . . . . . . . . . . . . . . . . . . . . . . 5.9 Fourier coefficients of u at t = 0.15 for an initial condition of the form (5.7). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Computed numerical energy for quintic NLS in Fig. 5.4 with N = 2 8

  Solution to the focusing DS II equation (3.2) for an initial condition of the form (5.8) with B = 0.1 for t = -3 and t = 0 in the first row and t = 3 and t = 6 below. . . . . . . . . . . . . . . . . . . . . . . . . 5.18 Evolution of max(|u| 2 ) and of the energy in dependence of time for an initial condition of the form (5.8) with B = 0.1. . . . . . . . . . . 5.19 Fourier coefficients of u at t = 6 for an initial condition of the form (5.8) with B = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.20 Evolution of max(|u| 2 ) in dependence of time for an initial condition of the form (5.8), B = 0.1, 0.5 . . . . . . . . . . . . . . . . . . . . . . 5.21 Contour plot for a solution to the focusing DS II equation (3.2) for

5 . 25

 525 Evolution of max(|u| 2 ) and the numerically computed energy Efor a solution to the focusing DS II equation (3.2) for an initial condition of the form (5.9) with κ = 1.1. . . . . . . . . . . . . . . . . . . . . . . 5.26 Fourier coefficients of the solution to the focusing DS II equation (3.2) for an initial condition of the form (5.9) with κ = 1.1 at t = 0. . . . . 5.27 Solution to the focusing DS II equation (3.2) for an initial condition of the form (5.10) with C = 1.1 for t = 0.075 and t = 0.15 in the first row and t = 0.225 and t = 0.3 below. . . . . . . . . . . . . . . . . . . 5.28 Evolution of max(|u| 2 ) and the numerically computed energy for an initial condition of the form (5.10) with C = 1.1. . . . . . . . . . . . . 5.29 Fourier coefficients of solution to the focusing DS II equation (3.2) for an initial condition of the form (5.10) with C = 1.1 at t = 0.15. . . 5.30 Evolution of max(|u| 2 ) in dependence of time, for an initial condition of the form (5.10) with C = 0.9. . . . . . . . . . . . . . . . . . . . . . 5.31 Solution to the focusing DS II equation (3.2) for an initial condition of the form (5.10) with C = 0.9, N t = 2000 for t = 0.075 and t = 0.15 in the first row and t = 0.225 and t = 0.3 below. . . . . . . . . . . . . 5.32 Fourier coefficients of the solution to the focusing DS II equation (3.2) at t = 0.15 for an initial condition of the form (5.10) with C = 0.9. . . 5.33 Solution to the focusing DS II equation (3.2) for an initial condition of the form (5.11) with D = 0.1 for t = 0.075 and t = 0.15 in the first row and t = 0.225 and t = 0.3 below . . . . . . . . . . . . . . . . . . . 5.34 Evolution of max(|u| 2 ) and the numerically computed energy in dependence of time for the solution to the focusing DS II equation (3.2) for an initial condition of the form (5.11) with D = 0.1. . . . . . . . . 5.35 Fourier coefficients of the solution to the focusing DS II equation (3.2) at t = 0.15 for an initial condition of the form (5.11) with D = 0.1. . 5.36 Evolution of max(|u| 2 ) and the numerically computed energy for the solution to the focusing DS II equation (3.2) for an initial condition of the form (5.11) with D = 0.5. . . . . . . . . . . . . . . . . . . . . . 5.37 Evolution of max(|u| 2 ) and the numerically computed energy E in dependence of time for a solution to the focusing DS II equation (3.2)

6. 1

 1 Behavior of |u| 2 for t = 0.065 and t = 0.13 and t = 0.2484 for the Ozawa solution and time evolution of max(|u num | 2 ) and of u numu exact 2 , DCRK method. . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Behavior of the Fourier Coefficients of u at t = 0.13 for the Ozawa solution, DCRK method. . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Evolution of max(|u| 2 ) and of log10(∆E) in dependence of time, for

Figure 1 . 1 : 1 Figure 1 . 2 :

 11112 Figure 1.1: ω in dependence of ξ for the 1-d wave equation, with α = 1

Figure 1 . 3 :

 13 Figure 1.3: Evolution of the solution to the linearized gKdV equation, with u 0 = -sech 2 (x).

Figure 1 . 4 :

 14 Figure 1.4: Time evolution of the solution of (1.25) with u(0, x) =sech 2 x and n = 1.

  [START_REF] Davey | On three-dimensional packets of surface waves[END_REF] with u(0, x) = -sech 2 (x) and ν = 0.1 in Fig.1.5.

Figure 1 . 5 :

 15 Figure 1.5: Time evolution of the solution of (1.29) with u(0, x) = -sech 2 (x) and ν = 0.1.

Figure 1 . 6 :

 16 Figure 1.6: Dispersive shock waves for KdV in small dispersion limit.

Figure 1 . 7 :

 17 Figure1.7: KdV and asymptotic solution with u 0 = -sech 2 (x) and = 0.01, see[START_REF] Grava | Numerical solution of the small dispersion limit of Korteweg de Vries and Whitham equations[END_REF].

Figure 1 . 8 :

 18 Figure 1.8: Solution to the defocusing cubic NLS equation in semiclassical limit with u(0, x) = e -x 2 , = 0.5 and t max = 1.

Figure 1 . 10 :

 110 Figure 1.10: Time evolution of the solution of (1.32) with u(0, x) = sech (x) and = 0.1, blow-up occurs.

Figure 1 . 11 :

 111 Figure 1.11: Zaitsev solution to the KP I equation, α = 1 and β = 0.5, at t = 0.

Figure 1 . 12 :

 112 Figure 1.12: Doubly periodic hyperelliptic KP II solution of genus 2 with parameters: b = 1, λ = 0.15, bλ 2 + d = -1, µ 1 = µ 2 = 0.25, ν 1 = -ν 2 = 0.25269207053125, ω 1 = ω 2 = -1.5429032317052, and ϕ 1,0 = ϕ 2,0 = 0, at t = 0.

Figure 2 . 1 :

 21 Figure 2.1: Example of Aliasing effect (x for f 1 , o for f 2 ) on the discrete grid.

Figure 3 . 1 :

 31 Figure 3.1: Normalized L 2 -norm of the numerical error in the solution to the KP I equation with initial data given by the Zaitsev solution for several numerical methods, as a function of N t (top) and as a function of CPU time (bottom).

Figure 3 . 2 :

 32 Figure 3.2: Solution to the KP I equation for the initial data u 0 = -∂ x sech 2 (R) where R = x 2 + y 2 for several values of t.The computation is carried out with 2 11 × 2 9 points for (x, y) ∈ [-5π, 5π] × [-5π, 5π], = 0.1 and t ≤ 0.4. As a reference solution, we consider the solution calculated with the Hochbruck-Ostermann method with N t = 5000 time steps. The dependence of the normalized L 2 -norm of the difference between this reference solution and the numerical solution is shown in Fig.3.3 in dependence on the time step with a regression analysis and in dependence on the CPU time. Here we can see clearly the phenomenon of order reduction established analytically for parabolic systems by Hochbruck and Ostermann[START_REF] Hochbruck | Exponential Runge-Kutta Methods for semilinear parabolic Problems[END_REF]. In the stiff regime (here up to errors of order 10 -4 ) DCRK does not converge, the Integrating Factor method shows only first order behavior (as predicted in[START_REF] Hochbruck | Exponential Runge-Kutta Methods for semilinear parabolic Problems[END_REF]), and ETD methods

Figure 3 . 3 :

 33 Figure 3.3: Normalized L 2 -norm of the numerical error for the solution shown in Fig. 3.2 for several numerical methods, as a function of N t (top) and as a function of CPU time (bottom).

Figure 3 . 4 :

 34 Figure 3.4: Phenomenon of order reduction in Exponential Integrators for KP I in the small dispersion limit, ETD schemes of Fig. 3.3

Figure 3 . 5 :

 35 Figure 3.5: Normalized L 2 -norm of the numerical error for the time evolution of the doubly periodic solution to the KP II equation for several numerical methods, as a function of N t (top) and as a function of CPU time (bottom).

Figure 3 . 6 :

 36 Figure 3.6: Solution to the KP II equation for the initial data u 0 = -∂ x sech 2 (R) where R = x 2 + y 2 for several values of t.

Figure 3 . 7 :

 37 Figure 3.7: Normalized L 2 -norm of the numerical error for the solution in Fig. 3.6 for several numerical methods as a function of N t (top) and as a function of CPU time (bottom).

  log10(Nt)+b1 -a2*log10(Nt)+b2

Figure 3 . 8 :

 38 Figure 3.8: Phenomenon of order reduction in Exponential Integrators for KPII in the small dispersion limit (see Fig. 3.7).

Figure 3 . 10 :

 310 Figure 3.10: Normalized L 2 -norm of the numerical error for several numerical methods for the situation shown in Fig. 3.9 as a function of N t (top) and of CPU time (bottom).

Figure 3 . 11 :

 311 Figure 3.11: Solution to the focusing DS II equation for the initial data u 0 = exp(-R 2 ) where R = x 2 + ηy 2 , η = 0.1 and = 0.1 for several values of t.

Figure 3 . 12 :

 312 Figure 3.12: Normalized L 2 -norm of the numerical error for the example of Fig. 3.11 for several numerical methods as a function of N t (top) and of CPU time (bottom).

Figure 3 . 13 :

 313 Figure 3.13: Non-conservation of the numerically computed L 2 -norm of the solution to the problem considered in Fig. 3.3 in dependence on the time step (top) and the Fourier coefficients for the final time (bottom).

Figure 3 .

 3 Figure 3.14: Non-conservation of the numerically computed L 2 -norm of the solution to the problem considered in Fig. 3.7 in dependence on the time step (top) and the Fourier coefficients for the final time (bottom).

Figure 3 .

 3 Figure 3.15: Non-conservation of the numerically computed L 2 -norm of the solution to the problem considered in Fig. 3.10 in dependence on the time step (top) and the Fourier coefficients for the final time (bottom).

Figure 3 .

 3 Figure 3.16: Non-conservation of the numerically computed L 2 -norm of the solution to the problem considered in Fig. 3.12 in dependence on the time step (top) and the Fourier coefficients for the final time (bottom).

Figure 4 . 1 :

 41 Figure 4.1: Serial Computation

Figure 4 . 2 :

 42 Figure 4.2: Example of three different parallel architectures to work on 4 different processes, P 0 , P 1 , P 2 and P 3 ; the double arrow representing a network.

Figure 4 . 4 :

 44 Figure 4.4: Example of data distribution for N x = N y = 8, n p = 4

Figure 4 . 5 :

 45 Figure 4.5: Reconstructed transposed matrix

Figure 4 . 6 :

 46 Figure 4.6: Scaling of the code

Figure 5 . 2 :

 52 Figure 5.2: Fourier coefficients for the situation in Fig. 5.1 at t = 6.

Figure 5 . 3 :

 53 Figure 5.3: Solution to the focusing quintic NLS (5.4) for the initial data u 0 = 1.8i exp(-x 2 ) with N = 2 12 on the left and N = 2 15 on the right for t > t c .

Figure 5 . 4 :

 54 Figure 5.4: Numerically computed energy for the situation studied in Fig. 5.3 for N = 2 12 on the left and N = on the right for several values of N t . At the blowup, the energy jumps.

Figure 5 . 5 :

 55 Figure 5.5: Fourier coefficients for the solution in Fig. 5.3 close to the critical and at a later time for N = 2 on the left and N = 2 15 on the right for N t = 10 4 .

7 )Remark 5 . 2

 752 (a = 1 and b = -4 in (1.2)). As for the quintic NLS in 1 + 1 dimensions, we always trace the conserved energy for DS II (3.2).The computation is carried out with N x = N y = 215 , L x = L y = 20, and N t = 1000 respectively N t = 4000; we show the solution at different times in Fig.5.6. The difference of the Ozawa mass and the computed L 2 norm on the periodic setting is of the order of 9 * 10 -5 . The time evolution of max x,y |u(x, y, t)| 2 and the difference between the numerical and the exact solution can be seen in Fig.5.7 (the critical time t c is not on the shown grid, thus the solution is always finite on the grid points). The code continues to run after the critical time, but the numerical solution obviously no longer represents the Ozawa solution. The numerically computed energy jumps at the blow up time as can be seen in Fig.5.8. The Fourier coefficients at t = 0.15 are shown in Fig.5.9. Despite the Gibbs phenomenon the Fourier coefficients for the initial data decrease to 10 -8 . Spatial resolution is still satisfactory at half the blowup time. The jump of the computed energy at blowup is dependent on sufficient spatial resolution as can be seen in Fig.5.10 for the example of the quintic NLS of Fig.5.4 and the Ozawa solution in Fig.5.8. For low resolution blow-up can be still clearly recognized from the computed energy, but the energy does not stay on the level at blow-up.

Figure 5 . 6 :

 56 Figure 5.6: Solution to the focusing DS II equation (3.2) for t = 0.075 and t = 0.15 in the first row and t = 0.225 and t = 0.3 below for an initial condition of the form (5.7).

Figure 5 . 7 :Figure 5 . 8 :

 5758 Figure 5.7: Time evolution of max(|u num | 2 ) and of u num -u exact 2 for the situation in Fig. 5.6.

Figure 5 . 9 :

 59 Figure 5.9: Fourier coefficients of u at t = 0.15 for an initial condition of the form (5.7).

Figure 5 . 10 :

 510 Figure 5.10: Computed numerical energy for quintic NLS in Fig. 5.4 with N = 2 8 and for the Ozawa solution in Fig. 5.8 with N x = N y = 2 12 .

  3) denoted by u l u(x, y, -6) = Au l , where A ∈ R is a scaling factor. The computations are carried out with N x = N y = 2 14 points for x × y ∈ [-50π, 50π] × [-50π, 50π] and t ∈ [-6, 6]. For A = 1.1, and N t = 1000, we observe a blowup of the solution at t c ∼ 1.6. The time evolution of max x,y |u(x, y, t)| 2 and of the energy is shown in Fig. 5.11. The maximum of |u| 2 in Fig.5.11 is clearly smaller than in the case of the Ozawa solution. This is due to the lower resolution in time which is used for this computation. Nevertheless, the jump in the energy is obviously present. The Fourier coefficients at t = 0 can be seen in Fig.5.12.

Figure 5 . 11 :

 511 Figure 5.11: Evolution of max(|u| 2 ) and the numerically computed energy in dependence of time for a solution to the focusing DS II equation (3.2) for an initial condition of the form u(x, y, -6) = 1.1u l .

Figure 5 . 12 :

 512 Figure 5.12: Fourier coefficients at t = 0 for a solution to the focusing DS II equation (3.2) for an initial condition of the form u(x, y, -6) = 1.1u l .

Figure 5 . 13 :

 513 Figure 5.13: Fourier coefficients at t = 6 for a solution to the focusing DS II equation (3.2) for an initial condition of the form u(x, y, -6) = 1.1u l .

  5.20.

Figure 5 . 14 :

 514 Figure 5.14: Evolution of max(|u| 2 ) and the numerically computed energy in dependence of time for a solution to the focusing DS II equation (3.2) for an initial condition of the form u(x, y, -6) = 0.9u l .

Figure 5 . 15 :

 515 Figure 5.15: Solution to the focusing DS II equation (3.2) for an initial condition of the form u(x, y, -6) = 0.9u l for t = -3 and t = 0 in the first row and t = 3 and t = 6 below.

Figure 5 . 16 :

 516 Figure 5.16: The Fourier coefficients at t = 0 of the solution to the focusing DS II equation (3.2) for an initial condition of the form u(x, y, -6) = 0.9u l .

Figure 5 . 17 :

 517 Figure 5.17: Solution to the focusing DS II equation (3.2) for an initial condition of the form (5.8) with B = 0.1 for t = -3 and t = 0 in the first row and t = 3 and t = 6 below.

Figure 5 . 18 :

 518 Figure 5.18: Evolution of max(|u| 2 ) and of the energy in dependence of time for an initial condition of the form (5.8) with B = 0.1.

Figure 5 . 19 :

 519 Figure 5.19: Fourier coefficients of u at t = 6 for an initial condition of the form (5.8) with B = 0.1.

5 Figure 5 . 20 :

 5520 Figure 5.20: Evolution of max(|u| 2 ) in dependence of time for an initial condition of the form (5.8), B = 0.1, 0.5

Figure 5 . 21 :

 521 Figure 5.21: Contour plot for a solution to the focusing DS II equation (3.2) for an initial condition of the form (5.9) with κ = 0.9 for different times.

Figure 5 . 22 :

 522 Figure 5.22: Evolution of max(|u| 2 ) and the numerically computed energy in dependence of time for the focusing DS II equation (3.2) for an initial condition of the form (5.9) with κ = 0.9.

-Figure 5 . 23 :Figure 5 . 24 :

 523524 Figure 5.23: Fourier coefficients of the solution to the focusing DS II equation (3.2) for an initial condition of the form (5.9) with κ = 0.9 at t = 0.

Figure 5 . 25 :

 525 Figure 5.25: Evolution of max(|u| 2 ) and the numerically computed energy Efor a solution to the focusing DS II equation (3.2) for an initial condition of the form (5.9) with κ = 1.1.

Figure 5 . 26 :

 526 Figure 5.26: Fourier coefficients of the solution to the focusing DS II equation (3.2) for an initial condition of the form (5.9) with κ = 1.1 at t = 0.

Figure 5 . 27 :

 527 Figure 5.27: Solution to the focusing DS II equation (3.2) for an initial condition of the form (5.10) with C = 1.1 for t = 0.075 and t = 0.15 in the first row and t = 0.225 and t = 0.3 below.

Figure 5 . 28 :

 528 Figure 5.28: Evolution of max(|u| 2 ) and the numerically computed energy for an initial condition of the form (5.10) with C = 1.1.

Figure 5 . 29 :

 529 Figure 5.29: Fourier coefficients of solution to the focusing DS II equation (3.2) for an initial condition of the form (5.10) with C = 1.1 at t = 0.15.

Figure 5 . 30 :

 530 Figure 5.30: Evolution of max(|u| 2 ) in dependence of time, for an initial condition of the form (5.10) with C = 0.9.

Figure 5 . 31 :

 531 Figure 5.31: Solution to the focusing DS II equation (3.2) for an initial condition of the form (5.10) with C = 0.9, N t = 2000 for t = 0.075 and t = 0.15 in the first row and t = 0.225 and t = 0.3 below.

Figure 5 . 32 :

 532 Figure 5.32: Fourier coefficients of the solution to the focusing DS II equation (3.2) at t = 0.15 for an initial condition of the form (5.10) with C = 0.9.

Figure 5 . 33 :Figure 5 . 34 :

 533534 Figure 5.33: Solution to the focusing DS II equation (3.2) for an initial condition of the form (5.11) with D = 0.1 for t = 0.075 and t = 0.15 in the first row and t = 0.225 and t = 0.3 below . energy indicating blowup at the time t c ∼ 0.2332. The Fourier coefficients at t c = 0.15 in Fig. 5.35 show that the wanted spatial resolution is achieved.

Figure 5 . 35 :

 535 Figure 5.35: Fourier coefficients of the solution to the focusing DS II equation (3.2) at t = 0.15 for an initial condition of the form (5.11) with D = 0.1.

  u(x, y, 0) = 2 exp -i(x 2 -(νy) 2 ) 1 + x 2 + (νy) 2 , (5.12) i.e., a deformation in the y-direction. The computations are carried out with N x = N y = 2 15 points for x × y ∈ [-20π, 20π] × [-20π, 20π] and t ∈ [0, 0.3].

Figure 5 . 36 :

 536 Figure 5.36: Evolution of max(|u| 2 ) and the numerically computed energy for the solution to the focusing DS II equation (3.2) for an initial condition of the form (5.11) with D = 0.5.

Figure 5 . 37 :

 537 Figure 5.37: Evolution of max(|u| 2 ) and the numerically computed energy E in dependence of time for a solution to the focusing DS II equation (3.2) for an initial condition of the form (5.12) with ν = 0.9.

Figure 5 . 38 :

 538 Figure 5.38: Fourier coefficients of the solution to the focusing DS II equation (3.2) for an initial condition of the form (5.12) with ν = 0.9 at t = 0.

Figure 5 . 39 :

 539 Figure 5.39: Evolution of max(|u| 2 ) and the numerically computed energy E for a solution to the focusing DS II equation (3.2) for an initial condition of the form (5.12) with ν = 1.1.

Figure 5 . 40 :

 540 Figure 5.40: Fourier coefficients of the solution to the focusing DS II equation (3.2) for an initial condition of the form (5.12) with ν = 1.1 at t = 0.15.

Figure 6 . 1 :

 61 Figure 6.1: Behavior of |u| 2 for t = 0.065 and t = 0.13 and t = 0.2484 for the Ozawa solution and time evolution of max(|u num | 2 ) and of u num -u exact 2 , DCRK method.

Figure 6 . 2 :

 62 Figure 6.2: Behavior of the Fourier Coefficients of u at t = 0.13 for the Ozawa solution, DCRK method.

6 . 3 .Figure 6 . 5 :

 6365 Figure 6.5: Time evolution of max(|u| 2 ) and of the computed energy, for an initial condition of the form u(x, y, -3) = 0.9u l , DCRK method.

Figure 6 . 6 :

 66 Figure 6.6: Behavior of the Fourier Coefficients of u at t = 6 for an initial condition of the form u(x, y, -3) = 0.9u l , DCRK method.

Figure 6 . 7 :

 67 Figure 6.7: Time evolution of the maximum of |u| 2 for an initial condition of the form (6.1), for A = 0.1, 0.3, 0.5.

Figure 6 . 8 :

 68 Figure 6.8: Behavior of |u| 2 at different times for initial data of the form (6.2) with η = 0.1.

Figure 6 . 9 :

 69 Figure 6.9: Time evolution of max(|u| 2 ) and of the computed energy for an initial data of the form (6.2) with η = 0.1.

6 . 19 .

 619 The comparison of the evolution of max(|u| 2 ) for an initial data of the form (6.2) for several values of η is also shown in Fig.6.19.

Figure 6 . 10 :

 610 Figure 6.10: Behavior of the Fourier Coefficients of u at t max = 0.6 for initial data of the form (6.2) with η = 0.1.

Figure 6 . 11 :

 611 Figure 6.11: Behavior of |u| 2 at different times for initial data of the form (6.2) with η = 0.3.

Figure 6 . 12 :

 612 Figure 6.12: Time evolution of max(|u| 2 ) and of the computed energy for initial data of the form (6.2) with η = 0.3.

Figure 6 . 13 :

 613 Figure 6.13: Behavior of the Fourier Coefficients of u at t max = 0.6 for initial data of the form (6.2) with η = 0.3.

Figure 6 . 14 :

 614 Figure 6.14: Behavior of |u| 2 at different times for initial data of the form (6.2) with η = 0.5.

Figure 6 . 15 :

 615 Figure 6.15: Time evolution of max(|u| 2 ) and of the computed energy for initial data of the form (6.2) with η = 0.5 (top) and η = 0.8 (bottom).

Figure 6 . 16 :

 616 Figure 6.16: Behavior of |u| 2 at different times for initial data of the form (6.2) with η = 0.8.

Figure 6 . 17 :

 617 Figure 6.17: Behavior of the Fourier Coefficients of u at t max = 0.6 for initial data of the form (6.2) with η = 0.5 (top) and η = 0.8 (bottom).

Figure 6 . 18 :

 618 Figure 6.18: Contours for different η for initial data of the form (6.2).

2 (

 2 log(|û(k, t)|) = A(t) -α(t) log(k) -δ(t)k (6.15)The error err on the fit interval k min < k < k max satisfieserr 2 = k min <k i <kmax (log(|û(k i , t)|) -(A(t) -α(t) log(k i ) -δ(t)k i )) linear in the fit parameters.

  3, 0.5. . . . . . . . . . . . . . . . . . . . . . 6.8 Behavior of |u| 2 at different times for initial data of the form (6.2) with η = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.9 Time evolution of max(|u| 2 ) and of the computed energy for an initial data of the form (6.2) with η = 0.1. . . . . . . . . . . . . . . . . . . . 6.10 Behavior of the Fourier Coefficients of u at t max = 0.6 for initial data of the form (6.2) with η = 0.1. . . . . . . . . . . . . . . . . . . . . . . 6.11 Behavior of |u| 2 at different times for initial data of the form (6.2) with η = 0.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.12 Time evolution of max(|u| 2 ) and of the computed energy for initial data of the form (6.2) with η = 0.3. . . . . . . . . . . . . . . . . . . . 6.13 Behavior of the Fourier Coefficients of u at t max = 0.6 for initial data of the form (6.2) with η = 0.3. . . . . . . . . . . . . . . . . . . . . . . 6.14 Behavior of |u| 2 at different times for initial data of the form (6.2) with η = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.15 Time evolution of max(|u| 2 ) and of the computed energy for initial data of the form (6.2) with η = 0.5 (top) and η = 0.8 (bottom). . . . 6.16 Behavior of |u| 2 at different times for initial data of the form (6.2) with η = 0.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.17 Behavior of the Fourier Coefficients of u at t max = 0.6 for initial data of the form (6.2) with η = 0.5 (top) and η = 0.8 (bottom). . . . . . . 6.18 Contours for different η for initial data of the form (6.2). . . . . . . . 6.19 Time evolution of max(|u| 2 ) and of the computed energy for initial Butcher table of a RK scheme with s stages . . . . . . . . . . . . . . 2.2 Butcher tables of (from left to right): Gauss method (of order 4 with s = 2), Radau I and Radau II method (of order 3 with s = 2). . . . . 2.3 Butcher table of a s-stage IF Runge-Kutta method, and of IFRK4 . . 2.4 Fourth order method of Cox-Matthews . . . . . . . . . . . . . . . . . 2.5 Fourth order method of Krogstad . . . . . . . . . . . . . . . . . . . . 2.6 Fourth order method of Hochbruck and Ostermann . . . . . . . . . .

	List of Tables
	2.1 12

data of the form (6.2) with η = 1 (top) and with different values of η (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 2

 2

.1 is not singular and if a sj = b j , for j = 1..s, and a i1 = b 1 , for i = 1..s,

(2.35)

then the method is L-stable. The method is stiffly stable if only the second condition is satisfied.

Thus, the Gauss and Radau IRK methods are preferred over other IRK schemes to solve stiff problems. Butcher tables of the Gauss method (of order 4 with s = 2), Radau I and II method (of order 3 with s = 2) are given in Table

2

.2, and we have the following result, see

[START_REF] Hairer | Solving ordinary differential equations[END_REF]

:

Theorem 2.5 1.

The s-stage Gauss method is A-stable and of order 2s.

2. The s-stage

Radau I and Radau II methods are A-stable and of order (2s -1).

Table 2 .

 2 3: Butcher table of a s-stage IF Runge-Kutta method, and of IFRK4

In quantum mechanics the Schrödinger equation plays the same fundamental role as Newton's equations of motion in classical mechanics and the Maxwell equations in the classical theory of electromagnetism.

In the framework of Hamiltonian mechanics, for solving equation (1.3) there should be n independent integrals of motion in involution.

Weak solutions being not unique, the notion of entropy solution was defined as follows: We call entropy for the Hopf equation (1.27) all pairs (U, F ) of functions C 1 from R to R such that U is a strictly convex function, and that F (u) = U (u)f (u), ∀ u ∈ R. A weak solution u to the Cauchy problem (1.27) is an entropy solution if for all entropies (U, F ) of the Hopf equation, it satisfies the entropy condition ∂ ∂t U (u) + ∂ ∂x F (u) ≤ 0 in the sense of distributions.

This was the first case for which the Riemann-Hilbert problem formulation of the IST scheme was inadequate.

All results in this section can be easily extended to an interval of different length than

2π.

Very smooth functions have very rapidly decaying Fourier coefficients, leading to the rapid convergence of the Fourier series, whereas discontinuous functions have very slowly decaying Fourier coefficients, causing the Fourier series to converge very slowly.

In more general cases, J admits a Jordan canonical form, so that (2.19) is obtained from (2.20) with M being the Jordan canonical form of J.

In our case this happens after the semi-discretization in space of dispersive PDEs.

The exact solutions described in Chap. 1 are obtained for = 1.

This work was granted access to the HPC resources of CCRT/IDRIS under the allocation

2011- [x2011106628] made by GENCI (Grand Equipement National de Calcul Intensif).

Acknowledgements

Tout d'abord, je tiens à remercier une des personnes sans qui je n'aurai jamais pensé faire une thèse, le Professeur Vladimir Matveev. Son enthousiasme et sa soif de connaissance communicatifs sont un bel exemple qui mérite le respect.

Je remercie également mon directeur de thèse, le Professeur Christian Klein, d'une part pour avoir accepté de travailler avec moi, et surtout pour m'avoir donné l'occasion de découvrir le monde de la recherche, avec ses bons et ses mauvais aspects, ainsi que l'opportunité de voyager lors de collaborations scientifiques, comme en Nouvelle Zélande, ou aux Etats Unis. Je n'oublierai jamais cette experience unique de ma vie, toutes ces découvertes tellement enrichissantes. J'en profite pour remercier M. Patrick Dutartre, sans qui je n'aurai pas encore réalisé tout ce que j'ai appris durant cette thèse, et a quel point j'ai changé, grandi, muri.

Au cours de ces trois ans, j'ai eu l'occasion de rencontrer nombre de personnes que je voudrai remercier pour diverses raisons. Tout d'abord, le Professeur Joerg Frauendiener, qui a eu la gentillesse de m'acceuillir en Nouvelle Zélande, et ce, à deux reprises; l'ambiance amicale et chaleureuse qui règne dans son équipe est une bouffée d'air frais dont on devrait s'inspirer en France. Un petit coucou à Florian et Lenette, merci de votre acceuil et de votre gentillesse. Parmi les nombreux thésards que j'ai rencontré, je voudrai remercier en particulier Georges, Phil, Nicolas, Katharina, Affaf, pour toutes les discussions intéressantes que nous avons eu; l'échange de nos expériences respectives m'a beaucoup apporté. Je remercie aussi Philippe, Muriel et Caroline de l'IMB, pour nos discussions, ainsi que Gauthier et Ahmed, pour nos meeting d'anglais.

Je remercie également M. Benson Muite, pour m'avoir invitée a une formation sur le calcul parallèle aux Etats-Unis, ainsi qu'à plusieurs conférences. Je tiens a remercier aussi M. Hervé Leferrant, pour le soutien et l'intérêt qu'il a porté à mon travail durant ces trois ans, ainsi que M. Olivier Politano pour toutes les formations en informatique qu'il a dispensé a l'universite, et sa disponibilité. Enfin, je remercie les personnes qui ont accepté de faire partie de mon jury de thèse, le Pr. Mariana Haragus, le Pr. Paul Matthews, le Pr. Vladimir Matveev, et tout

Stiff problems are those for which explicit methods are inefficient. and we will see why in the next section.

Absolute Stability and Stiff Problems

Consider for instance the solution to the very simple problem, y = -µy, with µ > 0, let y(t) = y(0)e -µt . If we now discretize it by using the simplest explicit numerical method, the explicit Euler scheme, we have

We thus have y n = (1-µ∆ t ) n y 0 . So, in particular, if ∆ t > 2 µ , (y n ) n∈N becomes unbounded, whereas the exact solution tends to 0 when n → ∞. In such cases, the study of the long time behavior of the solution is crucial to obtain satisfying approximations. We thus study in this context the absolute stability. For this study, we assume ∆ t fixed, and following Dahlquist (1963), we consider the model problem y = λy = f (t, y), λ ∈ C.

(2.22)

For the continuous problem, it is clear that y(t) tends to 0 when t → ∞ for all initial conditions, as soon as Re λ < 0. A numerical method is said to be A-stable if it satisfies the same property:

Definition 2.1 A-Stability

A method is said to be A-stable if, for any initial data y 0 , the numerical solution y n of (2.22) with Re λ < 0, tends to 0 as n → ∞ for any step size ∆ t .

Numerical schemes applied to the test case (2.22) are written as

where R : C → C is called the stability function of the method.

Remark 2.1

For homogeneity reasons, it is natural that parameters λ and ∆ t appear in the numerical schemes only in the product form λ∆ t .

Setting z = λ∆ t , clearly a numerical method will be A-stable if Re λ < 0 ⇒ |R(z)| < 1, which gives the following definition

Definition 2.2 A-Stability region

The set straight lines with slope a = 3.95 for the Integrating Factor method, a = 4.01 for DCRK, a = 3.99 for Krogstad's ETD scheme, a = 3.99 for the Cox-Matthews scheme, a = 3.95 for the Hochbruck-Ostermann scheme, and a = 3.86 for the time splitting method. The second order splitting scheme shows the expected convergence rate and performs very well for lower precision. For smaller time steps, the advantage of the fourth order schemes is more pronounced. Apparently the system is 'stiff' for the IRK4 scheme since it only shows second order behavior. Notice that the time splitting scheme reaches its maximal precision around 10 -8 , a behavior which was already noticed in [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger Equation[END_REF] for the study of the nonlinear Schrödinger equation in the small dispersion limit. It appears that this behavior is due to resonances of errors of the split equations, but the identification of the precise reason will be the subject of further research. The same effect is observed for second order splitting for smaller time steps than shown in Fig. 3.10. We conclude that the ETD schemes perform best in this context.

We choose N x = N y = 2 14 and L x = L y = 50, with ξ = 0, η = -1, z 0 = 1 and c = 1. The large values for L x and L y are necessary to ensure that the solution is small at the boundaries of the computational domain to reduce Gibbs phenomena. The difference for the mass of the lump and the computed mass on this periodic setting is of the order of 6 * 10 -5 . The initial data for t = -6 are propagated with N t = 1000 time steps until t = 6. In Fig. 5.1 contours of the solution at different times are shown. Here and in the following we always show closeups of the solution. The actual computation is done on the stated much larger domain. In this paper we will always show the square of the modulus of the complex solution for ease of presentation. The time dependence of the L 2 norm of the difference between the numerical and the exact solution can be also seen there. The numerical error is here mainly due to the lack of resolution in time. Since the increase in the number of time steps is computationally expensive, a fourth order scheme is very useful in this context. The spatial resolution can be seen from the modulus of the Fourier coefficients at the final time of computation t = 6 in Fig. 5.2. It decreases to 10 -6 , thus essentially the value for the initial data. For computational speed considerations we always use double precision which because of finite precision arithmetic give us a range of 15 orders of magnitude. Since function values computed using the split step method were for most of the computation of order 1, and less than 5,000, rounding errors allow for a precision of 10 -14 when less than 2 15 × 2 15 Fourier modes are used. When more modes than 2 15 × 2 15 were used, we found a reduction in precision. Despite the algebraic falloff of the solution we have a satisfactory spatial resolution because of the large computational domain and the high resolution. The modulational instability does not show up in this and later examples before blowup.

machine precision, despite the modulational instability. For A = 0.9, we recover exactly the behavior observed with the splitting method ( the initial pulse travels in the same direction as the exact solution, but loses speed and height and is broadened), see Fig. 6.5 and the computed energy is obviously smoother. The Fourier coefficients are shown in Fig. 6.6 at t = 6, they again reach machine precision.

We compare the time evolution of the maximum of |u| 2 from initial data of the form

for different values of A in Fig. 6.7 (perturbation of the Ozawa solution with a Gaussian). The corresponding times of blow up are 0.2335 for A = 0.1, 0.1994 for A = 0.3 and 0.1664 for A = 0.5. It thus appears proportional to the perturbation.

In view of the clear better performance of the DCRK method, compared to the splitting

Conclusion/Outlook

We mainly discussed in this work the efficient numerical integration of dispersive nonlinear equations, where the stiffness of the systems is related to the linear part. In particular we studied two (2+1)-dimensional equations, the Kadomtsev-Petviashvili equation and the Davey-Stewartson system for which so far no analytical results exist for cases of particular interest, as the behavior of the solutions in the small dispersion limit, or the integrability of the corresponding dispersionless equations, and also blow-up phenomena.

Where theories fail, numerical simulations become the only way to investigate these problems. They allow to validate existing theories, can help in he formulation of conjectures and provide directions of future research by advancing. This applies of course only to carefully done numerical simulations with control of potential errors as in this thesis.

We began by comparing several time stepping schemes to identify the most efficient method for the different studied PDEs. We found that it is not possible to select a single optimal time integration scheme for all situations. In fact, we were forced to identify a scheme which performs best for one particular equation and even for a particular regime, mainly in terms of stiff and non-stiff regimes.

From convergence tests and the (non)-conservation of quantities such as the mass and energy for the equations under consideration, we found that fourth order time stepping schemes can be efficiently used to deal with stiff systems and dispersive shocks. In particular, Driscoll's composite RK method is generally very efficient if the studied system is not too stiff, but fails to converge for strong stiffness. We also observed that the order reduction phenomenon, first studied by Hochbruck and Ostermann for parabolic PDEs, can be considerable for IF schemes in the stiff regime, but less so for ETD schemes in the case of the here considered hyperbolic equations. The Hochbruck-Ostermann method which performs in general best, requires in practice more CPU time and suffers, as others methods from a small order reduction in stiff regimes.

The behavior of solutions of such equations in the small dispersion limit (or semiclassical limit) was the main topic of this work. An understanding of the semiclassical limit of a 2+1-dimensional equations is even more challenging than the 1+1-dimensional cases for which a complete understanding for general analytic initial data is sometimes still missing.

A first question to investigate would be the existence of soliton solutions for these systems.

In [START_REF] Ghidaglia | Nonexistence of Travelling Wave Solutions to Nonelliptic Nonlinear Schrödinger Equations[END_REF], it was shown that the hyperbolic-elliptic DS systems does not support travelling wave solutions except for a specific range of the parameters α , β, that includes the focusing DS IIcase. However no conjecture is given for the case β < 0 , α ∈ (0, -β) and α + β/2 = 0, which remains an open problem. A numerical study of this could be very useful.

We could also give indications of the behavior of the solutions of the defocusing DS equations in the non-integrable cases, since the existence of solutions for defocusing DS equations is only known for the integrable version, see [START_REF] Fokas | On the solvability of the N-wave, the Davey-Stewartson and the Kadomtsev-Petviashmli equation[END_REF].

Thirdly, the behavior of the time-splitting scheme has to be understood theoretically.

To this end, some experiments have to be done, to find in which cases these methods can be efficiently used, and how they can be improved.

Last, but not least, the study of the dispersionless DS system (6.5

) not yet achieved, is the focus of future work. A powerful tool to determine the appearence of singularities could be the analyticity strip method, introduced by Sulem and Sulem and Frisch in [START_REF] Sulem | Tracing complex singularities with spectral methods[END_REF]. The basic idea is to obtain the width of the analyticity strip, from the asymptotic behavior of the Fourier transform û(k, t) (in one spatial dimension), or from the angle averaged energy spectrum defined by On the other hand, consider an analytic function u(z) in one spatial dimension defined on R or on a periodic domain with singularities at the complex locations z j , in the neighborhood of which it behaves as u(z) ∼ (z -z j ) µ j . The behavior of the Fourier transform of such a function for k → +∞ is dominated by the singularity in the upper half-space closest to the real axis that is not a multiple pole. Denoting this singularity of exponent µ by z * = x * + iδ, one has û(k) ∼ |k| -(µ+1) e -kδ e ix * k , k → +∞ (6.11)

A derivation of this property can be found for instance in [START_REF] Carrier | Functions of a complex variable: Theory and technique[END_REF]. It requires that u(z) does not grow faster than an exponential as k → +∞ and that the singularities are isolated. Equation (6.11) implies in particular that for a function of one variable, the width of the