Numerical Study of Blowup in the Davey-Stewartson System

NUMERICAL STUD Y OF BLO WUP IN THE D A VEY-STEW AR TSON SYSTEM C. KLEIN, B. MUITE, AND K. ROIDOT Abstract. Nonlinear dispersive partial differential equations suc h as the non- linear Schr¨ odinger equations can hav e solutions that blo w-up. W e n umerically study the long time b ehavior and potential blo wup of solutions to the focusing Dav ey-Stew artson I I equation b y analyzing p erturbations of the lump and the Ozaw a exact solutions. It is sho wn in this w a y that the lump is unstable to both blowup and dispersion, and that blowup in the Oza w a solution is generic. 1. Introduction The Da vey-Stew artson (DS) system models the evolution of w eakly nonlinear w ater w a ves that trav el predominan tly in one direction for whic h the w av e ampli- tude is mo dulated slowly in tw o horizon tal directions [8], [9]. It is also used in plasma ph ysics [20, 21], to describe the evolution of a plasma under the action of a magnetic field. The DS system can b e written in the form (1) i∂ t u + ∂ xx u − α∂ y y u + 2 ρ  Φ + | u | 2  u = 0 , ∂ xx Φ + β ∂ y y Φ + 2 ∂ xx | u | 2 = 0 , where α , β and ρ tak e the v alues ± 1, and where Φ is a mean field. The DS equations can be seen as a tw o-dimensional nonlinear Schr¨ odinger (NLS) equation with a nonlo cal term if the equation for Φ can b e solv ed for given b oundary conditions. They are classified [12] according to the ellipticity or hyperb olicity of the operators in the first and second line in (1). The case α = β is completely in tegrable [1] and th us provides a 2 + 1-dimensional generalization of the integrable NLS equation in 1 + 1 dimensions with a cubic nonlinearity . The integrable cases are elliptic- h yp erbolic called DS I, and the h yp erbolic-elliptic called DS I I. F or b oth there is a fo cusing ( ρ = − 1) and a defo cusing ( ρ = 1) version. The complete integrabilit y of the DS equation implies that it has an infinite num b er of conserved quantities, for Key words and phr ases. Dav ey-Stew artson systems, split step, blow-up. W e thank J.-C. Saut for helpful discussions and hints. CK and KR were supp orted by the pro ject F roM-PDE funded by the European Researc h Council through the Adv anced Inv estigator Grant Scheme, the Conseil R´ egional de Bourgogne via a F ABER grant, the Marie-Curie IRSES program RIMMP and the ANR via the program ANR-09-BLAN-0117-01. BM is grateful to the Institut de Math ´ ematiques de Bourgogne for their hospitality and financial supp ort by the CNRS where part of this work was completed. Part of the computational requirements of this research were supported b y an allocation of adv anced computing resources pro vided b y the National Science F oundation. This w ork w as gran ted access to the HPC resources of CCR T/IDRIS under the allocation 2011106628 made by GENCI, Kraken at the National Institute for Computational Sciences and Hopper at NERSC. Computational supp ort was also provided by CRI de Bourgogne and SCREMS NSF DMS-1026317. 1 2 C. KLEIN, B. MUITE, AND K. ROIDOT instance the L 2 norm and the energy E [ u ( t )] := 1 2 Z T 2  | ∂ x u ( t, x, y ) | 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 . DS reduces to the cubic NLS in one dimension if the p oten tial is indep enden t of y , and if Φ satisfies certain b oundary conditions (for instance rapidly decreasing at infinit y or p eriodic). In the following, w e will only consider the case DS I I ( α = 1) since the mean field Φ is then obtained by inv erting an elliptic op erator. The non- in tegrable elliptic-elliptic DS is very similar to the 2 + 1 dimensional NLS equation, see for instance [12] and [6] for numerical sim ulations, and is therefore not studied here. There exist many explicit solutions for the integrable cases which thus allow to address the question ab out the long time b ehavior of solutions for given initial data. F or the famous Korteweg-de V ries (KdV) equation, it is known that general initial data are either radiated aw a y or asymptotically decomp ose into solitons. The DS I I system and the t w o-dimensional in tegrable generalization of KdV known as the Kadom tsev-P etviash vili I (KP I) equation ha v e so-called lump solutions, a t wo-dimensional soliton whic h is lo calized in all spatial directions with an algebraic falloff to wards infinity . F or KP I it was shown [2] that small initial data asymptot- ically decomp ose into radiation and lumps. It is conjectured that this is also true for general initial data. F or the defo cusing DS I I global existence in time was sho wn by F ok as and Sung [10] for solutions of certain classes of Cauch y problems. These initial data will simply disperse. The situation is more in volv ed for the focusing case. P elino vski and Sulem [24] show ed that the lump solution is sp ectrally unstable. In addition the fo cusing NLS equations in 2 + 1 dimensions with cubic nonlinearit y ha v e the critical dimension, i.e., solutions from smo oth initial data can ha v e blowup . This means that the solutions lose after finite time the regularity of the initial data, a norm of the solution or of one of its deriv atives b ecomes infinite. F or fo cusing NLS equations in 2 + 1 dimensions, it is known that blowup is p ossible if the energy of the initial data is greater than the energy of the ground state solution, see e.g. [27] and references therein, and [19] for an asymptotic description of the blowup profile. F or the fo cusing DS I I equation Sung [28] ga v e a smallness condition on the F ourier transform F [ u 0 ] of the initial data to establish global existence in time for solutions to Cauch y problems (2) ||F [ u 0 ] || L 1 ||F [ u 0 ] || L ∞ < π 3 2 √ 5 − 1 2 ! 2 ∼ 5 . 92 . . . . with initial data u 0 ∈ L p , 1 ≤ p < 2 with a F ourier transform F [ u 0 ] ∈ L 1 ∩ L ∞ . It is not kno wn 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 pap er are not in this class, w e cannot pro vide further insigh t in to this question. An explicit solution with blo wup for lump-like initial data w as given b y Ozaw a [22]. It has an L ∞ blo wup in one p oin t ( x c , y c , t c ) and is analytic for all other v alues of ( x, y , t ). It is unkno wn whether this is the typical blowup b eha vior for the fo cusing DS I I equation. NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 3 F rom the p oin t of view of applications, a blowup of a solution do es not mean that the studied equation is not relev ant in this con text. It just indicates the limit of the used approximation. It is th us of particular interest, not only in mathematics, but also in physics, since it shows the limits of the applicability of the studied mo del. This breakdown of the mo del will also in general indicate how to amend the used appro ximations. In view of the op en analytical questions concerning blowup in DS I I solutions, w e study the issue in the present pap er numerically , which is a highly non-trivial problem for several reasons: first DS is a purely dispersive equation whic h means that the introduction of n umerical dissipation has to b e av oided as m uc h as p ossible to preserve disp ersiv e effects suc h as rapid oscillations. This makes the use of sp ectral methods attractive since they are known for minimal n umerical dissipation and for their excellent approximation prop erties for smo oth functions. But the algebraic falloff of b oth the lump and the Ozaw a solution leads to strong Gibbs phenomena at the b oundaries of the computational domain if the solutions are p eriodically contin ued there. W e will nonetheless use F ourier spectral methods b ecause they also allow for efficient time in tegration algorithms whic h should b e ideally of high order to av oid a p ollution of the F ourier co efficien ts due to n umerical errors in the time integration. An additional problem is the mo dulational instability of the focusing DS II equa- tion, i.e., a self-induced amplitude mo dulation of a contin uous wa ve propagating in a nonlinear medium, with subsequent generation of lo calized structures, see for instance [3, 7, 11] for the NLS equation. Th us to address n umerically questions of stability and blowup of its solutions, high resolution is needed which cannot b e ac hieved on single pro cessor computers. Therefore we use parallel computers to study the related questions. The use of F ourier sp ectral metho d is also very con- v enient in this context, since for a parallel sp ectral code only existing optimized serial FFT algorithms are necessary . In addition such codes are not memory in- tensiv e, in con trast to other approaches such as finite difference or finite elemen t metho ds. The first numerical studies of DS were done by White and W eideman [32] using F ourier sp ectral metho ds for the spatial co ordinates and a second order time splitting sc heme. Besse, Mauser and Stimming [6] used essen tially a paral- lel v ersion of this co de to study the Ozaw a solution and blowup in the focusing elliptic-elliptic DS equation. McConnell, F ok as and Pelloni [18] used W eideman’s co de to study n umerically DS I and DS I I, but did not hav e enough resolution to get conclusiv e results for the blowup in p erturbations of the lump in the fo cusing DS I I case. In this pap er we rep eat some of their computations with considerably higher resolution. W e use a parallelized version of a fourth order time splitting scheme which was studied for DS in [17]. Obviously it is non-trivial to decide numerically whether a solution blo ws up or whether it just has a strong maxim um. T o allo w to mak e nonetheless reliable statements, we p erform a series of tests for the numerics. First w e test the co de on known exact solutions with algebraic falloff, the lump and the Oza wa solution. W e establish that energy conserv ation can b e used to judge the qualit y of the numerics if a sufficient spatial resolution is givem. It is shown that the splitting co de contin ues to run in general beyond a p oten tial blowup whic h mak es it difficult to decide whether there is blowup. W e argue at examples for the quintic NLS in 1 + 1 dimensions (which is kno wn to hav e blowup solutions) 4 C. KLEIN, B. MUITE, AND K. ROIDOT and the Oza wa solution that energy conserv ation is a reliable indicator in this case since the energy of the solution changes completely after a blowup, whereas it will b e in accordance with the numerical accuracy after a strong maxim um. Th us w e repro duce well known blo wup cases in this wa y and establish with the energy conserv ation a criterion to ensure the accuracy of the numerics also in unknown cases. Then we study p erturbations of the lump and the Ozaw a solution to see when blowup is actually observed. The pap er is organized as follows: in section 2, w e describ e the numerical co de and its parallelization, and study as an example the lump solution. In section 3 w e numerically study blowup in the fo cusing 1+1-dimensional quintic NLS and the Oza wa solution for DS I I. In section 4 we discuss p erturbations of the lump, and in section 5 perturbations of the Oza wa solution. In section 6 we give some concluding remarks. 2. Numerical methods In this pap er we are in terested in the n umerical study of solutions to the focusing DS I I equation for initial data with algebraic falloff to w ards infinity in all spatial directions. This algebraic decrease of the initial data and consequently of the solution for all times is in principle not an ideal setting for the use of F ourier metho ds since the perio dic contin uation of the function at the b oundaries of the computational domain will lead to Gibbs phenomena. Nonetheless there are several reasons for the use of F ourier metho ds in this case: First the DS equation is a purely dispersive PDE, and we are in terested in disp ersiv e effects. Thus it imp ortan t to use n umerical metho ds that introduce as little nu- merical dissipation as p ossible, and sp ectral metho ds are esp ecially effective in this con text. F urthermore the discrete F ourier transform can b e efficien tly computed with a fast F ourier tr ansform (FFT). In addition F ourier metho ds allo w to use splitting techniques for the tim e integration as explained below in a very efficient w ay . Last but not least the fo cusing DS I I equation is kno wn to hav e a modulation instabilit y which makes the use of high resolution necessary , esp ecially close to the blo wup situations w e wan t to study . This instability leads to an artificial increase of the high wa ve num b ers which even tually breaks the co de, if not enough spatial resolution is pro vided (see for instance [16] for the focusing NLS equation). It is not p ossible to reac h the necessary resolution on single pro cessors which mak es a paral- lelization of the co de obligatory . As explained b elo w, this can b e conv eniently done for 2-dimensional (even for 3-dimensional) F ourier transformations where the task of the 1-dimensional FFTs is performed simultaneously by sev eral pro cessors. This reduces also the memory requirements p er pro cessor ov er alternative approac hes. 2.1. Splitting Metho ds. Splitting metho ds are v ery efficient if an equation can b e split into tw o or more equations which can b e directly integrated. They are un- conditionally stable. The motiv ation for these metho ds is the T rotter-Kato formula [31, 15] (3) lim n →∞  e − tA/n e − tB /n  n = e − t ( A + B ) where A and B are certain unbounded linear operators, for details see [15]. In particular this includes the cases studied by Bagrinovskii and Godunov in [5] and b y NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 5 Strang [26]. F or hyperb olic equations, first references are T app ert [30] and Hardin and T app ert [14] who introduced the split step metho d for the NLS equation. The idea of these methods for an equation of the form u t = ( 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 num b ers that represent frac- tional time steps. Y oshida [33] ga v e an approac h which produces split step metho ds of any even order. The DS equation can b e split into i∂ t u = ( − ∂ xx u + α∂ y y u ) , ∂ xx Φ + α∂ y y Φ + 2 ∂ xx  | u | 2  = 0 , (4) i∂ t u = − 2 ρ  Φ + | u | 2  u, (5) whic h are explicitly integrable, the first tw o in F ourier space, equation (5) in ph ys- ical space since | u | 2 and thus Φ is constant in time for this equation. Con vergence of second order splitting along these lines was studied in [6]. W e use here fourth order splitting as giv en in [33] and already studied in [17] for the DS II equation. In the latter reference, it w as shown that this scheme is v ery efficient in this context. The metho d is conv enien t for parallel computing, b ecause of easy co ding (lo ops) and low memory requirements. Notice that the splitting metho d in the form (5) conserv es the L 2 norm: the first equation implies that its solution in F ourier space is just the initial condition (from the last time step) multiplied b y a factor e iφ with φ ∈ R . Thus the L 2 norm is constan t for solutions to this equation b ecause of P arsev al’s theorem. The second equation as men tioned conserves the L 2 norm exactly . Th us the used splitting sc heme has the conserv ation of the L 2 norm implemented. As we will show in the follo wing, this does not guarantee the accuracy of the n umerical solution since other conserv ed quantities as the energy the conserv ation of which is not implemented migh t not b e numerically conserved. In fact we will sho w that the numerically computed energy provides a v alid indicator of the quality of the numerics. 2.2. P arallelization of the co de. Since high resolution is needed to numerically examine the fo cusing DS I I equation, the co de is parallelized to reduce the wall clo c k time required to run the sim ulation and to allo w the problem to fit in memory . The runs typically used N x = N y = 2 15 , where N x and N y denote the num b er of F ourier mo des in x and y resp ectiv ely . The parallelization is done by a slab domain decomp osition. The grid p oin ts are given by x n = 2 π nL x N x , y m = 2 π mL y N y , so that the numerical solution is in the computational domain x × y ∈ [ − L x π , L x π ] × [ − L y π , L y π ] . In the computations, L x = L y is c hosen large enough so that the n umerical solution is small at the boundaries, and hence a n umerical solution on a perio dic domain can b e considered as a go od approximation to the solution on an un b ounded do- main. The appro ximate solution u is represen ted by an N x × N y matrix, which is distributed among the MPI pro cesses (note that each MPI pro cess uses a single core). F or programming ease and for the efficiency of the F ourier transform, N x and N y are chosen to be p o wers of t w o. The n umber of MPI pro cesses, n p is chosen 6 C. KLEIN, B. MUITE, AND K. ROIDOT to divide N x and N y p erfectly , so that each pro cess holds N x × N y /n p elemen ts of u , for example pro cess i holds the elements u (1 : N x , ( i − 1) N y n p + 1 : i N y n p ) in the global arra y . T o a v oid p erforming global F ourier transforms whic h are in- efficien t, the arra y is transp osed once all the one dimensional F ourier transforms in the x direction hav e b een done. Since the data is evenly distributed among the MPI pro cesses, this transp ose is efficiently implemented using MPI ALLTOALL [13]. After the transp ose, the F ourier transform ˆ u is distributed on the pro cesses so that pro cess i holds the elements corresp onding to ˆ u (( i − 1) N x n p + 1 : iN x n p , 1 : N y ) , on whic h the second set of one dimensional FFTs can be done. The one dimensional FFTs were done using FFTW 3.0, FFTW 3.1 and FFTW 3.2 1 whic h are close to optimal on x86 architectures and allow the resulting program to b e portable but still simple. 2.3. Lump solution of the fo cusing DS I I equation. T o test the p erformance 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 fo cusing DS I I equation has solitonic solutions which are regular for all x, y , t , and which are lo calized with an algebraic falloff tow ards infinity , known as lumps [4]. The single lump is given b y (6) u ( x, y , t ) = 2 c exp  − 2 i ( ξ x − η y + 2( ξ 2 − η 2 ) t )  | x + 4 ξ t + i ( y + 4 η t ) + z 0 | 2 + | c | 2 where ( c, z 0 ) ∈ C 2 and ( ξ , η ) ∈ R 2 are constants. The lump mov es with constan t v elo cit y ( − 2 ξ , − 2 η ) and decays as ( x 2 + y 2 ) − 1 for x, y → ∞ . W e choose N x = N y = 2 14 and L x = L y = 50, with ξ = 0 , η = − 1 , z 0 = 1 and c = 1. The large v alues for L x and L y are necessary to ensure that the solution is small at the b oundaries of the computational domain to reduce Gibbs phenomena. The difference for the mass of the lump and the computed mass on this p eriodic setting is of the order of 6 ∗ 10 − 5 . The initial data for t = − 6 are propagated with N t = 1000 time steps un til t = 6. In Fig. 1 con tours of the solution at differen t times are sho wn. Here and in the follo wing w e alw ays sho w closeups of the solution. The actual computation is done on the stated muc h larger domain. In this pap er w e will alwa ys show the square of the mo dulus of the complex solution for ease of presen tation. The time dep endence of the L 2 norm of the difference b et ween the n umerical and the exact solution can b e also seen there. The numerical error is here mainly due to the lack of resolution in time. Since the increase in the n umber of time steps is computationally exp ensiv e, a fourth order sc heme is very useful in this context. The spatial resolution can b e seen from the mo dulus of the F ourier co efficien ts at the final time of computation t = 6 in Fig. 2. It decreases to 10 − 6 , th us essen tially the v alue for the initial data. F or computational sp eed c onsiderations w e alw a ys use double precision whic h because of finite precision arithmetic give us a range of 15 orders of magnitude. Since function 1 http://www.fftw.org NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 7 Figure 1. Con tours of | u | 2 on the left and a plot of || u exact − u num || 2 on the right in dep endence of time for the solution to the fo cusing DS I I equation (1) for lump initial data (6). v alues computed using the split ste p metho d were for most of the computation of order 1, and less than 5,000, rounding errors allo w for a precision of 10 − 14 when less than 2 15 × 2 15 F ourier mo des are used. When more modes than 2 15 × 2 15 w ere used, we found a reduction in precision. Despite the algebraic falloff of the solution w e ha ve a satisfactory spatial resolution b ecause of the large computational domain and the high resolution. The mo dulational instability do es not show up in this and later examples b efore blowup. Figure 2. F ourier co efficien ts for the situation in Fig. 1 at t = 6. 3. Blowup for the quintic NLS in 1 + 1 dimensions and the focusing DS I I It is known that fo cusing NLS equations can hav e solutions with blowup, if the nonlinearit y exceeds a critical v alue depending on the spatial dimension. F or the 1 + 1 dimensional case, the critical nonlinearit y is quin tic, for the 2 + 1 dimensional it is cubic, see for instance [27] and references therein. Thus the fo cusing DS II equation can ha v e solutions with blowup. In this section w e will first study numerically blo wup for the 1 + 1 dimensional quintic NLS equation, and then numerically evolv e initial data for a known exact blo wup solution to the fo cusing DS I I equation due to Ozaw a [22]. W e discuss some p eculiarities of the fourth order splitting scheme in this context. 8 C. KLEIN, B. MUITE, AND K. ROIDOT 3.1. Blo wup for the quintic one-dimensional NLS. The fo cusing quin tic NLS in 1 + 1 dimensions has the form (7) i∂ t u + ∂ xx u + | u | 4 u = 0 , where u ∈ C dep ends on x and t (w e consider again solutions p erio dic in x ). This equation is not completely integrable, but assuming the solution is in L 2 , has conserv ed L 2 norm and, provided the solution u ∈ H 2 , a conserved energy , (8) E [ u ] = Z R  1 2 | ∂ x u | 2 − 1 6 | u | 6  dx. It is known that initial data with negative energy blo w up for this equation in finite time, and that the b ehavior close to blo wup is giv en in terms of a solution to an ODE, see [19]. As discussed in sect. 2.1, the splitting scheme we are using here has the prop ert y that the L 2 norm is conserv ed. Thus the qualit y of the numerical conserv ation of the L 2 norm giv es no indication on the accuracy of the n umerical solution. Ho wev er as discussed in [16], conserv ation of the n umerically computed energy giv es a v alid criterion for the qualit y of the numerics: in general it ov erestimates the L ∞ n umerical error by tw o orders of magnitude at the typically aimed at precisions. If we consider as in [25] for the quintic NLS the initial data u 0 ( x ) = 1 . 8 i exp( − x 2 ), the energy is negative. W e compute the solution with L x = 5 and N x = 2 15 with N t = 10 4 time steps. The result can b e seen in Fig. 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 low er spatial resolution N = 2 12 ). The initial data ï 2 ï 1 0 1 2 0 0.1 0.2 0 100 200 300 x t |u| 2 Figure 3. Solution to the fo cusing quintic NLS (7) for the initial data u 0 = 1 . 8 i exp( − x 2 ) with N = 2 12 on the left and N = 2 15 on the right for t > t c . clearly get fo cused to a strong maximum, but the co de does not break. W e note that this is in contrast to other fourth order schemes tested for 1 + 1 dimensional NLS equations in [16], whic h t ypically pro duce an o verflo w close to the blo wup. But clearly the solution shows spurious oscillations after the time t c ∼ 0 . 155. In fact the numerically computed energy , whic h will alwa ys b e time-dependent due to unav oidable n umerical errors, will be completely changed after this time. W e consider (9) ∆ E =     1 − E ( t ) E (0)     , NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 9 where E ( t ) is the n umerically computed energy (8) and get for the example in Fig. 3 the b eha vior shown in Fig. 4. A t the presumed blowup at t c ∼ 0 . 155 as in [25], the energy jumps to a completely different v alue. Thus this jump can and will b e used to indicate blowup. T o illustrate the effects of a low er resolution in time and space imp osed by hardware limitations for the DS computations, we show this quantit y for several resolutions in Fig. 4. If a low er resolution in time is used as in some of the DS examples in this pap er, the jump is slightly smo othed out. But the plateau is still reached at essentially the same time which indicates blowup. Th us a lack of resolution in time in the given limits will not b e an obstacle to identify a p ossible singularit y . The reason for this is the use of a fourth order scheme that allows to tak e larger time steps. W e will present computations with different resolutions to illustrate the steep ening of the energy jump as ab o ve if this is within the limitations imp osed b y the hardware. 0 0.05 0.1 0.15 0.2 ï 12 ï 10 ï 8 ï 6 ï 4 ï 2 0 2 4 t log 10 6 E N t = 1000 N t = 3000 N t = 10000 0 0.05 0.1 0.15 0.2 ï 15 ï 10 ï 5 0 5 t log 10 6 E N t = 1000 N t = 3000 N t = 10000 Figure 4. Numerically computed energy for the situation studied in Fig. 3 for N = 2 12 on the left and N = 2 15 on the righ t for sev eral v alues of N t . At the blo wup, the energy jumps. W e show the mo dulus of the F ourier co efficien ts for N = 2 12 and N = 2 15 b efore and after the critical time in Fig. 5. It can b e seen that the solution is well resolv ed b efore blo wup in the latter case, and that the singularity leads to oscillations in the F ourier coefficients. A lack of spatial resolution as for N = 2 12 in Fig. 5 triggers the mo dulation instabilit y close to the blowup and at later times as can b e seen from the mo dulus of the F ourier co efficien ts that increase for larger wa ven um b ers. Therefore we alwa ys aim at a sufficient resolution in space even for times close to a blo wup. After this time the mo dulation instabilit y will b e present in the spurious solution pro duced by the splitting scheme as we will show for an example. Remark 3.1. Stinis [25] has r e c ently c ompute d singular solutions to the fo cusing quintic nonline ar Schr¨ odinger e quation in 1 + 1 dimensions. This e quation has solutions in L ∞ L 2 that may not b e unique for given smo oth initial data and that may exhibit blowup of the L ∞ H 1 norm. F ol lowing T ao [29] , Stinis [25] has use d a sele ction criteria to pick a solution after the blow up time of the L ∞ H 1 norm. They suggest that ‘mass’ is eje cte d (which me ans that the L 2 norm is change d) at times wher e the L ∞ H 1 norm blows up. The splitting scheme studie d her e in c ontr ast pr o duc es a we ak solution with a differ ent ener gy sinc e the L 2 norm c onservation is built in. 10 C. KLEIN, B. MUITE, AND K. ROIDOT ï 500 0 500 ï 5 0 5 t=0.134 log 10 |v| ï 500 0 500 ï 5 0 5 t=0.136 log 10 |v| ï 500 0 500 ï 5 0 5 k t=0.2 log 10 |v| ï 4000 ï 2000 0 2000 4000 ï 20 ï 10 0 10 t=0.134 log 10 |v| ï 4000 ï 2000 0 2000 4000 ï 20 ï 10 0 10 log 10 |v| t=0.136 ï 4000 ï 2000 0 2000 4000 ï 20 ï 10 0 10 log 10 |v| k t=0.2 Figure 5. F ourier co efficien ts for the solution in Fig. 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 . 3.2. Blo wup in the Oza w a solution. F or the fo cusing DS I I equation, an exact solution was given by Ozaw a [22] which is in L 2 for all times with an L ∞ blo wup in finite time. W e can summarize his results as follows: Theorem 3.1 (Ozaw a) . L et ab < 0 and T = − a/b . Denote by u ( x, y , t ) the function define d by (10) u ( x, y , t ) = exp  i b 4( a + bt ) ( x 2 − y 2 )  v ( X, Y ) a + bt wher e (11) v ( X, Y ) = 2 1 + X 2 + Y 2 , X = x a + bt , Y = y a + bt Then, u is a solution of (1) with (12) k u ( x, y , t ) k 2 = k v ( X, Y ) k 2 = 2 √ π and (13) | u ( t ) | 2 → 2 π δ when t → T . wher e δ is the Dir ac me asur e. W e thus consider initial data of the form (14) u ( x, y , 0) = 2 exp  − i ( x 2 − y 2 )  1 + x 2 + y 2 ( a = 1 and b = − 4 in (10)). As for the quintic NLS in 1 + 1 dimensions, w e alwa ys trace the conserved energy for DS I I (1). The computation is carried out with N x = N y = 2 15 , L x = L y = 20, and N t = 1000 resp ectiv ely N t = 4000; we show the solution at different times in Fig. 6. The difference of the Oza w a mass and the computed L 2 norm on the p eriodic setting is of the order of 9 ∗ 10 − 5 . The time evolution of max x,y | u ( x, y , t ) | 2 and the difference b etw een the n umerical and the exact solution can be seen in Fig. 7 (the critical time t c is not on the sho wn grid, thus the solution is alw a ys finite on the grid p oints). The co de contin ues to run after the critical time, but the numerical solution ob viously no longer represents the Ozaw a solution. The numerically computed energy jumps at the blow up time NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 11 Figure 6. Solution to the fo cusing DS I I equation (1) for t = 0 . 075 and t = 0 . 15 in the first row and t = 0 . 225 and t = 0 . 3 b elow for an initial condition of the form (14). Figure 7. Time ev olution of max ( | u num | 2 ) and of k u num − u exact k 2 for the situation in Fig. 6. as can b e seen in Fig. 8. The F ourier co efficien ts at t = 0 . 15 are shown in Fig. 9. Despite the Gibbs phenomenon the F ourier co efficien ts for the initial data decrease to 10 − 8 . Spatial resolution is still satisfactory at half the blo wup time. Remark 3.2. The jump of the c ompute d ener gy at blowup is dep endent on sufficient sp atial r esolution as c an b e se en in Fig. 10 for the example of the quintic NLS of Fig. 4 and the Ozawa solution in Fig. 8. F or low r esolution blow-up c an b e stil l cle arly r e c o gnize d fr om the c ompute d ener gy, but the ener gy do es not stay on the level at blow-up. 12 C. KLEIN, B. MUITE, AND K. ROIDOT Figure 8. Numerically computed energy E ( t ) and ∆ E = | 1 − E ( t ) /E (0) | (9) for the situation in Fig. 6. Figure 9. F ourier co efficien ts of u at t = 0 . 15 for an initial con- dition of the form (14). 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 ï 12 ï 11 ï 10 ï 9 ï 8 ï 7 ï 6 ï 5 ï 4 ï 3 t log 10 6 E Figure 10. Computed numerical energy for quintic NLS in Fig. 4 with N = 2 8 and for the Ozaw a solution in Fig. 8 with N x = N y = 2 12 . 4. Per turba tions of the lump solution In this section w e consider p erturbations of the lump solution (6). First we propagate initial data obtained from the lump after multipli cation with some scalar factor. Then we consider a p erturbation with a Gaussian and a deformed lump. NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 13 4.1. P erturbation of the lump by rescaled initial data. W e first consider rescaled initial data from the lump (6) 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 p oin ts for x × y ∈ [ − 50 π , 50 π ] × [ − 50 π , 50 π ] and t ∈ [ − 6 , 6]. F or 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 sho wn in Fig. 11. The maxim um of | u | 2 in Fig. 11 is clearly smaller than in the case of the Oza w a solution. This is due to the lo w er resolution in time which is used for this computation. Nev- ertheless, the jump in the energy is ob viously presen t. The F ourier co efficien ts at Figure 11. Ev olution of max ( | u | 2 ) and the numerically computed energy in dep endence of time for a solution to the fo cusing DS I I equation (1) for an initial condition of the form u ( x, y , − 6) = 1 . 1 u l . t = 0 can b e seen in Fig. 12. They again decrease by almost 6 orders of magnitude. T o illustrate the mo dulational instabilit y at a concrete example, we show the F ourier co efficients after the critical time in Fig. 13. It can b e seen that the mo d- ulus of the co efficien ts of the high wa ven um bers increases instead of decreasing as to b e exp ected for smo oth functions. This indicates once more that the computed solution after the blowup time has to b e taken with a grain of salt. F or A = 0 . 9, the initial pulse trav els in the same direction as the exact solution, but loses sp eed and height and is broadened, see Fig. 14. It appears that this mo d- ified lump just disp erses asymptotically . The solution can be seen in Fig. 15. Its F ourier co efficien ts in Fig. 16 show that the resolution of the initial data is almost main tained. 4.2. P erturbation of the lump with a Gaussian. W e consider an initial con- dition of the form (15) u ( x, y , − 6) = u l + B exp( − ( x 2 + y 2 )) , B ∈ R . F or B = 0 . 1 and N t = 1000, we show the solution at different times in Fig. 17. The solution trav els at the same sp eed as b efore, but its amplitude v aries, growing and decreasing successively , see Fig. 18. The time evolution of the energy can b e seen in Fig. 18. There is no indication of blo wup in this example. The solution app ears 14 C. KLEIN, B. MUITE, AND K. ROIDOT Figure 12. F ourier co efficien ts at t = 0 for a solution to the fo cusing DS I I equation (1) for an initial condition of the form u ( x, y , − 6) = 1 . 1 u l . Figure 13. F ourier co efficien ts at t = 6 for a solution to the fo cusing DS I I equation (1) for an initial condition of the form u ( x, y , − 6) = 1 . 1 u l . Figure 14. Ev olution of max ( | u | 2 ) and the numerically computed energy in dep endence of time for a solution to the fo cusing DS I I equation (1) for an initial condition of the form u ( x, y , − 6) = 0 . 9 u l . to disp erse for t → ∞ . The F ourier co efficien ts at t = 6 in Fig. 19 show the wan ted spatial resolution. NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 15 Figure 15. Solution to the fo cusing DS I I equation (1) for an initial condition of the form u ( x, y , − 6) = 0 . 9 u l for t = − 3 and t = 0 in the first ro w and t = 3 and t = 6 b elo w. Figure 16. The F ourier co efficien ts at t = 0 of the solution to the fo cusing DS I I equation (1) for an initial condition of the form u ( x, y , − 6) = 0 . 9 u l . A similar b ehavior is observ ed if a larger v alue for the amplitude of the p ertur- bation is chosen, e.g., B = 0 . 5. 4.3. Deformation of the Lump. W e consider initial data of the form (16) u ( x, y , − 6) = u l ( x, κy , − 6) , i.e., a deformed (in y -direction) initial lump in this subsection. The computations are carried out with N x = N y = 2 14 p oin ts for x × y ∈ [ − 50 π , 50 π ] × [ − 50 π , 50 π ] and t ∈ [ − 6 , 6]. 16 C. KLEIN, B. MUITE, AND K. ROIDOT Figure 17. Solution to the fo cusing DS I I equation (1) for an initial condition of the form (15) with B = 0 . 1 for t = − 3 and t = 0 in the first ro w and t = 3 and t = 6 b elo w. Figure 18. Evolution of max ( | u | 2 ) and of the energy in dep en- dence of time for an initial condition of the form (15) with B = 0 . 1. F or κ = 0 . 9, the resulting solution loses sp eed and width as can b e seen in Fig. 20. Its heigh t and energy grow, but b oth stay finite, see Fig. 21. It is p ossible that the solution even tually blows up, but not on the time scales studied here. The F ourier co efficients at t = 0 in Fig. 22 show the wan ted spatial resolution. F or κ = 1 . 1, w e observe the opp osite b eha vior in Fig. 23. The solution tra v els with higher sp eed than the initial lump and is broadened. The energy do es not NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 17 Figure 19. F ourier co efficients of u at t = 6 for an initial condition of the form (15) with B = 0 . 1. −10 −5 0 5 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 Figure 20. Con tour plot for a solution to the fo cusing DS I I equa- tion (1) for an initial condition of the form (16) with κ = 0 . 9 for differen t times. Figure 21. Ev olution of max ( | u | 2 ) and the numerically computed energy in dep endence of time for the focusing DS I I equation (1) for an initial condition of the form (16) with κ = 0 . 9. sho w an y sudden c hange, see Fig. 24. It seems that the initial pulse will asymptot- ically disp erse. The F ourier co efficien ts at t = 0 in Fig. 25 show the wan ted spatial resolution. 18 C. KLEIN, B. MUITE, AND K. ROIDOT −200 −150 −100 −50 0 50 100 150 200 −8 −6 −4 −2 0 2 4 kx log 10 |v(kx,0)| Figure 22. F ourier co efficien ts of the solution to the fo cusing DS I I equation (1) for an initial condition of the form (16) with κ = 0 . 9 at t = 0. −15 −10 −5 0 5 10 15 −10 −8 −6 −4 −2 0 2 4 6 8 10 Figure 23. Con tour plot for a solution to the fo cusing DS I I equa- tion (1) for an initial condition of the form (16) with κ = 1 . 1 for differen t times. Figure 24. Ev olution of max ( | u | 2 ) and the numerically computed energy E for a solution to the fo cusing DS II equation (1) for an initial condition of the form (16) with κ = 1 . 1. NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 19 −200 −150 −100 −50 0 50 100 150 200 −8 −6 −4 −2 0 2 4 kx log 10 |v(kx,0)| Figure 25. F ourier co efficien ts of the solution to the fo cusing DS I I equation (1) for an initial condition of the form (16) with κ = 1 . 1 at t = 0. 5. Per turba tions of the Oza w a solution In this section we study as for the lump in the previous section v arious p ertur- bations of initial data for the Ozaw a solution to test whether blowup is generic for the fo cusing DS I I equation. 5.1. P erturbation of the Ozaw a solution b y m ultiplication with a scalar factor. W e consider initial data of the form (17) u ( x, y , 0) = 2 C exp  − i ( x 2 − y 2 )  1 + x 2 + y 2 , i.e., initial data of the Ozaw a solution multiplied by a scalar factor. The computa- tion is carried out with N x = N y = 2 15 p oin ts for x × y ∈ [ − 20 π , 20 π ] × [ − 20 π , 20 π ]. F or C = 1 . 1, and N t = 2000, w e show the behavior of | u | 2 at different times in Fig. 26. The time evolution of max x,y | u ( x, y , t ) | 2 and the numerically computed energy are shown in Fig. 27. W e observ e an L ∞ blo wup at the time t c ∼ 0 . 2210. The F ourier co efficien ts at t = 0 . 15 (b efore the blowup) in Fig. 28 show the wan ted spatial resolution. F or C = 0 . 9, the initial pulse gro ws until it reaches its maximal height at t = 0 . 2501, but there is no indication for blowup, see Fig. 29. The solution at different times can b e seen in Fig. 30. The F ourier co efficien ts in Fig. 31 show again that the wan ted spatial resolution is achiev ed. Th us for initial data given by the Ozaw a solution m ultiplied with a factor C , we find that for C > 1, blow up seems to o ccur b efore the critical time of the Ozaw a solution, and for C < 1 the solution grows until t = 0 . 25 but do es not blow up. Consequen tly the Ozaw a initial data seem to b e critical in this sense that data of this form with smaller norm do not blow up. 5.2. P erturbation of the Ozaw a solution with a Gaussian. W e consider an initial condition of the form (18) u ( x, y , 0) = 2 exp  − i ( x 2 − y 2 )  1 + x 2 + y 2 + D exp( − ( x 2 + y 2 )) . 20 C. KLEIN, B. MUITE, AND K. ROIDOT Figure 26. Solution to the fo cusing DS I I equation (1) for an initial condition of the form (17) with C = 1 . 1 for t = 0 . 075 and t = 0 . 15 in the first ro w and t = 0 . 225 and t = 0 . 3 b elo w. Figure 27. Ev olution of max ( | u | 2 ) and the numerically computed energy for an initial condition of the form (17) with C = 1 . 1. F or D = 0 . 1 and N t = 2000, we show the b eha vior of | u | 2 at differen t times in Fig. 32. The time evolution of max x,y | u ( x, y , t ) | 2 is sho wn in Fig. 33. W e observe a jump of the energy indicating blo wup at the time t c ∼ 0 . 2332. The F ourier co- efficien ts at t c = 0 . 15 in Fig. 34 sho w that the w anted spatial resolution is achiev ed. The same exp erimen t with D = 0 . 5 app ears again to sho w blow up, but at an earlier time t c ∼ 0 . 1659, see Fig. 35. Th us the energy added b y the p erturbation of the form D exp( − ( x 2 + y 2 )) seems NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 21 Figure 28. F ourier co efficien ts of solution to the fo cusing DS II equation (1) for an initial condition of the form (17) with C = 1 . 1 at t = 0 . 15. Figure 29. Ev olution of max ( | u | 2 ) in dep endence of time, for an initial condition of the form (17) with C = 0 . 9. to lead to a blowup before the critical time of the Ozaw a solution. This means that the blo wup in the Ozaw a solution is clearly a generic feature at least for initial data close to Ozaw a for the fo cusing DS I I equation. 5.3. Deformation of the Ozaw a solution. W e study deformations of Ozaw a initial data of the form (19) u ( x, y , 0) = 2 exp  − i ( x 2 − ( ν y ) 2 )  1 + x 2 + ( ν y ) 2 , i.e., a deformation in the y -direction. The computations are carried out with N x = N y = 2 15 p oin ts for x × y ∈ [ − 20 π , 20 π ] × [ − 20 π , 20 π ] and t ∈ [0 , 0 . 3]. F or ν = 0 . 9, we observe a maximum of the solution at t = 0 . 2441, see Fig. 36, follo wed b y a second maximum, but there is no indication of a blowup. Energy conserv ation is in principle high enough to indicate that the solution stays regular on the considered time scales. The F ourier coefficients at t = 0 . 15 in Fig. 37 show the w anted spatial resolution. 22 C. KLEIN, B. MUITE, AND K. ROIDOT Figure 30. Solution to the fo cusing DS I I equation (1) for an initial condition of the form (17) 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 b elo w. Figure 31. F ourier co efficien ts of the solution to the fo cusing DS I I equation (1) at t = 0 . 15 for an initial condition of the form (17) with C = 0 . 9. The situation is similar for ν = 1 . 1. The maxim um of the solution is observ ed at t = 0 . 2254, see Fig. 38, follow ed again b y a second maximum. Energy conserv ation app ears once more to rule out a blowup in this case. The F ourier co efficien ts at t = 0 . 15 in Fig. 39 again show the wan ted spatial resolution. NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 23 Figure 32. Solution to the fo cusing DS I I equation (1) for an initial condition of the form (18) with D = 0 . 1 for t = 0 . 075 and t = 0 . 15 in the first ro w and t = 0 . 225 and t = 0 . 3 b elo w . Figure 33. Ev olution of max ( | u | 2 ) and the numerically computed energy in dep endence of time for the solution to the focusing DS I I equation (1) for an initial condition of the form (18) with D = 0 . 1. 6. Conclusion In this pap er we ha ve numerically studied long time b eha vior and stability of exact solutions to the focusing DS I I equation with an algebraic falloff tow ards infinit y . W e hav e shown that the necessary resolution can be achiev ed with a parallelized version of a sp ectral co de. The spatial resolution as seen at the F ourier co efficien ts was alw a ys well b ey ond typical plotting accuracies of the order of 10 − 3 . F or the time integration we used an unconditionally stable fourth order splitting 24 C. KLEIN, B. MUITE, AND K. ROIDOT Figure 34. F ourier co efficien ts of the solution to the fo cusing DS I I equation (1) at t = 0 . 15 for an initial condition of the form (18) with D = 0 . 1. Figure 35. Ev olution of max ( | u | 2 ) and the numerically computed energy for the solution to the focusing DS I I equation (1) for an initial condition of the form (18) with D = 0 . 5. Figure 36. Ev olution of max ( | u | 2 ) and the numerically computed energy E in dep endence of time for a solution to the fo cusing DS I I equation (1) for an initial condition of the form (19) with ν = 0 . 9. sc heme. As argued in [16, 17], the numerically computed energy of the solution giv es a v alid indicator of the accuracy for sufficient spatial resolution. T o ensure the latter, we alwa ys presented the F ourier co efficien ts of the solution at a time b efore NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 25 Figure 37. F ourier co efficien ts of the solution to the fo cusing DS I I equation (1) for an initial condition of the form (19) with ν = 0 . 9 at t = 0. Figure 38. Ev olution of max ( | u | 2 ) and the numerically computed energy E for a solution to the fo cusing DS I I equation (1) for an initial condition of the form (19) with ν = 1 . 1. Figure 39. F ourier co efficien ts of the solution to the fo cusing DS I I equation (1) for an initial condition of the form (19) with ν = 1 . 1 at t = 0 . 15. a singularity app eared. In addition we sho w here that the numerically computed energy indicates blowup b y jumping to a different v alue in cases where the co de runs b ey ond a singularity in time. 26 C. KLEIN, B. MUITE, AND K. ROIDOT After testing the co de for exact solutions, the lump and the blo wup solution by Oza wa, w e sho w ed that both solutions are critical in the following sense: adding energy to it leads to a blowup for the lump, and an earlier blo wup time for the Ozaw a solution. F or initial data with less energy , no blowup was observed in b oth cases, the initial data asymptotically just seem to be disp ersed. This is in accordance with the conjecture in [18] that solutions to the fo cusing DS I I equations either blow up or disp erse. In particular the lump is unstable against b oth blowup and disp ersion, in con trast to the lump of the KP I equation that app ears to b e stable, see for instance [23]. Note that the p erturbations w e considered here test the nonlinear regime of the PDE for whic h so far no analytical results app ear to b e established. References [1] M. Ablowitz and R. Haberman , Nonline ar Evolution Equations in Two and Thre e Dimen- sions , Phys. Rev. Lett., 35 (1975), pp. 1185–8. [2] M. J. Ablowitz and A. Fokas , On the inverse sc attering and dir e ct line arizing tr ansforms for the Kadomtsev-Petviashvili e quation , Physics Letters A, 94 (1983), pp. 67–70. [3] G. Agra w al , Nonline ar Fiber Optics , Academic Press, San Diego, 2006. [4] V. Arkadiev, A. Pogrebkov, and M. Poliv anov , Inverse Sc attering Tr ansform Metho d and Soliton Solutions for Davey-Stewartson II Equation , Physica D: Nonlinear Phenomena, 36 (1989), pp. 189–197. [5] K. A. Bagrinovskii and S. Godunov , Differ enc e Schemes for multi-dimensional Pr oblems , Dokl.Acad. Nauk., 115 (1957), pp. 431–433. [6] C. Besse, N. Mauser, and H. Stimming , Numeric al Study of the Davey-Stewartson System , Math. Model. Numer. Anal., 38 (2004), pp. 1035–1054. [7] M. Cross and P. Hohenberg , Pattern formation outside of e quilibrium , Rev. Mo d. Phys., 65 (1993), p. 8511112. [8] A. Da vey and K. Stew ar tson , On Thr e e-dimensional Packets of Surfac e Waves , Pro c. R. Soc. Lond. A., 338 (1974), pp. 101–110. [9] V. Djordjevic and L. Redekopp , On Two-dimensional Packets of Capil larity-Gr avity Waves , J. Fluid Mech., 79 part 4 (1977), pp. 703–714. [10] A. Fokas and L. Sung , The Inverse sp e ctr al Metho d for the KP I Equation without the Zero Mass Constraint , Math. Proc. Camb. Phil. Soc., 125 (1999), pp. 113–138. [11] M. F orest and J. Lee , Ge ometry and mo dulation the ory for the p erio dic nonline ar Schr¨ odinger equation , in Oscillation Theory , Computation, and Metho ds of Comp ensated Compactness, Minneap olis, MN, 1985. The IMA V olumes in Mathematics and Its Applica- tions, vol. 2, Springer, New Y ork, 1986, pp. 35–69. [12] J.-M. Ghida glia and J.-C. Saut , On the initial Value Pr oblem for the Davey-Stewartson Systems , Nonlinearity , 3 (1990). [13] W. Gropp, E. Lusk, and A. Skjellum , Using MPI , MIT Press, Boston, 1999. [14] R. H. Hardin and F. D. T apper t , Applic ations of the Split-Step Fourier Metho d to the numeric al Solution of nonline ar and variable Coefficient Wave Equations , SIAM Rev., 15 (1973), p. 423. [15] T. Ka to , T r otter’s pr o duct formula for an arbitr ary p air of self-adjoint c ontr action semi- gr oups , vol. 3, Academic Press, Boston, 1978, pp. 185–195. [16] C. Klein , F ourth-Or der Time-Stepping for low Dispersion Kortewe g-de Vries and nonline ar Schr¨ odinger Equation , Electronic T ransactions on Numerical Analysis., 39 (2008), pp. 116– 135. [17] C. Klein and K. Roidot , F ourth or der time-stepping for Kadomtsev-Petviashvili and Davey- Stewartson equations , SIAM J. Sci. Comp., (2011). [18] M. McConnell, A. Fokas, and B. Pelloni , L o calise d c oher ent Solutions of the DSI and DSII Equations a numeric al Study , Mathematics and Computers in Simulation, 69 (2005), pp. 424–438. [19] F. Merle and P. Raphael , On universality of blow-up profile for L 2 critic al nonline ar Schr¨ odinger e quation , Inv entiones Mathematicae, 156 (2004), pp. 565–672. NUMERICAL STUDY OF BLOWUP IN THE DA VEY-STEW AR TSON SYSTEM 27 [20] K. Nishinari, K. Abe, and J. Sa tsuma , A new Typ e of Soliton Behaviour in a Two dimen- sional Plasma System , J. Phys. Soc. Jpn., 62 (1993), pp. 2021–2029. [21] , Multidimensional Behavior of an ele ctr ostaic Ion Wave in a magnetized Plasma , Phys. Plasmas, 1 (1994), pp. 2559–2565. [22] T. Oza w a , Exact Blow-Up Solutions to the Cauchy Pr oblem for the Davey-Stewartson Sys- tems , Proc. R. So c. Lond. A., 436 (1992), pp. 345–349. [23] D. Pelinovsky and C. Sulem , Eigenfunctions and Eigenvalues for a sc alar Riemann-Hilb ert Pr oblem asso ciate d to Inverse Sc attering , Comm un. Math. Phys., 208 (2000), pp. 713–760. [24] , Sp e ctr al dec omp osition for the dirac system asso ciate d to the dsii e quation , Inv. Prob., 16 (2000), pp. 59–74. [25] Stinis , Numeric al c omputation of solutions of the critic al nonline ar Schr¨ odinger e quation after the singularity , arXiv:1010.2246v1, (2010). [26] G. Strang , On the Construction and Comp arison of Differenc e Schemes , SIAM J. Numer. Anal., 5 (1968), pp. 506–517. [27] C. Sulem and P.-L. Sulem , The nonline ar Schr¨ odinger Equation , vol. 139 of Applied Math- ematical Sciences, Springer-V erlag, New Y ork, 1999. [28] L.-Y. Sung , L ong-Time Dec ay of the Solutions of the Davey-Stewartson II Equations , J. Nonlinear Sci., 5 (1995), pp. 433–452. [29] T. T ao , Global existenc e and uniqueness r esults for we ak solutions of the fo cusing mass- critic al nonline ar Schr¨ odinger e quation , Analysis and PDE, 2 (2009), pp. 61–81. [30] F. T apper t , Numeric al Solutions of the Kortewe g-de Vries Equation and its Gener alizations by the Split-Step Fourier Metho d , Lectures in Applied Mathematics, 15 (1974), pp. 215–216. [31] H. Trotter , On the Pr o duct of Semi-Gr oups of Op er ators , Proceedings of the American Mathematical Society , 10 (1959), pp. 545–551. [32] P. White and J. Weideman , Numerical Simulation of Solitons and Dromions in the Davey- Stewartson System , Math. Comput. Simul., 37 (1994), pp. 469–479. [33] H. Yoshida , Construction of higher Or der symple ctic Inte gr ators , Physics Letters A, 150 (1990), pp. 262–268. Institut de Ma th ´ ema tiques de Bourgogne, Universit ´ e de Bourgogne, 9 a venue Alain Sa v ar y, 21078 Dijon Cedex, France E-mail address : christian.klein@u-bourgogne.fr Dep ar tment of Ma thema tics, Uni versity of Michigan, 2074 East Hall, 530 Church Street, MI 48109, USA E-mail address : muite@umich.edu Institut de Ma th ´ ema tiques de Bourgogne, Universit ´ e de Bourgogne, 9 a venue Alain Sa v ar y, 21078 Dijon Cedex, France E-mail address : kristelle.roidot@u-bourgogne.fr