Numerical study of blow-up in solutions to generalized Kadomtsev-Petviashvili equations

We present a numerical study of solutions to the generalized Kadomtsev-Petviashvili equations with critical and supercritical nonlinearity for localized initial data with a single minimum and single maximum. In the cases with blow-up, we use a dynamic rescaling to identify the type of the singularity. We present a discussion of the observed blow-up scenarios.

1. Introduction. The celebrated Korteweg-de Vries (KdV) equation modelling wave phenomena in 1+1 dimensions in the limit of long wavelengths appears in many domains of application as hydrodynamics, acoustics, nonlinear optics, plasma physics, . . . . It has a 2 + 1 dimensional generalization for essentially one-dimensional wave phenomena with weak transverse effects, the Kadomtsev-Petviashvili (KP) equations [10] with a similar range of applications. Remarkably both KdV and KP are completely integrable. Therefore many explicit solutions as solitons are known for these equations. The Cauchy problem for the KP equations and the stability of certain exact solutions are actively studied, see for instance [16] for a recent review.
Since both KdV and KP appear in approximations for long wavelengths, the dispersion in these equations is too strong compared to what is observed in applications. A possible cure to this short coming is to shift the balance between nonlinearity and dispersion towards the nonlinearity. This leads for KdV to generalized KdV (gKdV) equations, and for KP to generalized KP (gKP) equations u t + u n u x + u xxx + λ∂ −1 x u yy = 0 (1.1) where n ∈ Q, n ≥ 1 and where λ = ±1; the antiderivative ∂ −1 x is defined via its Fourier symbol −i/k x . Note that gKP for n = 2 appears as a model for the evolution of sound waves in antiferromagnetic materials, see [6]. We will always consider asymptotically (for |x|, |y| → ∞) decreasing solutions in the following. The case λ = 1 is denoted gKP II, the one with λ = −1 gKP I. gKP I has a focusing effect, gKP II a defocusing one. If u is independent of y, the gKP equations reduce to the gKdV equations; for n = 1, one obtains the standard KP equations. The equations are only completely integrable in the latter case, but have conserved quantities as the mass, the L 2 norm of u, 2) and the energy Note that the equations (1.1) are invariant under rescalings of the form x → x/s, y → y/s 2 , t → t/s 3 and u → s 2/n u with s = const. Since the L 2 norm of u is invariant under this rescaling for n = 4/3, this case is also referred to as L 2 critical. For gKdV the critical exponent in this sense is n = 4.
The appearance of an antiderivative in (1.1) has the consequence that R ∂ yy u(x, y, t) dx = 0, ∀t > 0, (1.4) which can be easily seen by differentiating (1.1) with respect to x and integrating over the real line. This is even true if this constraint is not satisfied for the initial data u 0 (x, y). In [7,21] it was shown that the solution to a Cauchy problem not satisfying the constraint will not be smooth in time for t = 0. Numerical experiments in [17] indicate that the solution develops an infinite 'trench' the integral over which just ensures that (1.4) is fulfilled. This non-regularity of the solution for t → 0 is numerically a problem in the sense that it delimits the achievable accuracy.
To avoid such problems, we always consider here initial data that are x-derivatives of a rapidly decreasing function and thus satisfy (1.4); as a concrete example we will study in this paper the initial data u 0 (x, y) = β∂ xx exp(−(x 2 + y 2 )) (1.5) with β > 1.
Another consequence of the antiderivative in (1.1) is that solutions to Cauchy problems with initial data u 0 (x, y) in the Schwarz space S(R 2 ) of rapidly decreasing functions will not stay in this space unless u 0 (x, y) satisfies an infinite number of constraints. This can be already seen on the level of the linearized KP equation, see e.g. [2,17], where the Green's function implies a slow algebraic decrease in y for |y| → ∞. This leads to the formation of tails with an algebraic decrease towards infinity for generic localized initial data. The amplitude of these grows with time (see for instance [17]).
The gKP I equations are known to have localized travelling (in x-direction) wave solutions called lumps for n < 4. This corresponds to solutions to gKP of the form u(x, y, t) = Q(x − ct, y), c = const, where (the nontrivial) Q(z, y) satisfies − cQ zz + 1 n + 1 (Q n+1 ) zz + Q zzzz + λQ yy = 0. (1.6) The only explicit form of these solutions is known for KP I in terms of rational functions, see [23]. For general nonlinearities, solutions of (1.6) were studied in [4,5]. It was shown that they have an algebraic fall-off in x and y for |x|, |y| → ∞, i.e., that the solutions cannot decrease more rapidly than 1/(x 2 + y 2 ). The precise fall off rate was given in [8]. The lumps are unstable for n > 4/3, see [4,5]. There are no such solutions for gKP II. For KP I it was argued in [1] for small initial data that the solutions for large t decompose into lumps and radiation, a fact stressing the importance of lumps. It is conjectured that this behavior is also true for general initial data. Due to the instability of lumps for gKP I equations with n ≥ 4/3, an evolution of initial data into an array of lumps cannot be expected for this case for large t. But it is interesting to study which structures can be observed in the time evolution of localized initial data and which blow up. This is one of the goals of this paper. For gKdV, it was found numerically in [13] that initial data with sufficient mass can lead to the formation of several solitons, and that the one appearing first will finally blow up.
It is well known that there can be blow-up, i.e., a loss of regularity with respect to the initial data, in solutions to gKdV and gKP I equations with n ≥ 4 and n ≥ 4/3 respectively for certain initial data, see [19]. For the latter it was shown in [21] that the L 2 norm of u y blows up in these cases. It is unclear whether there can be blow-up in solutions to gKP II equations. A first numerical study of these issues has been presented in [16] 1 . In this paper we present a more detailed numerical analysis of these questions with the goal to identify the type of blow-up for localized initial data with a single minimum.
The paper is organized as follows: in Section 2 we present the used numerical methods and discuss a dynamic rescaling of the gKP equations. In section 3 we consider the L 2 critical case n = 4/3 and discuss various examples. In the cases with blow-up, we try to identify the type of the singularity.
The same analysis is performed in section 4 for the case n = 2 as an example for a supercritical situation. In section 5 we consider blow-up in gKP II solutions for n = 3 and n = 4. We add some concluding remarks in section 6.
2. Numerical methods. In this section we present the numerical methods to be used in this paper to integrate the gKP equations. The main tool will be a direct integration of the equations with a Fourier spectral method for the spatial coordinates, and a fourth order exponential time differencing (ETD) scheme for the time integration. We also discuss a dynamic rescaling of the equation in order to study blow-up cases in more detail.
2.1. Direct numerical integration. Fourier spectral approaches are generally the most efficient method for the numerical treatment of smooth periodic functions. The reason for this is the in practice exponential decrease of the Fourier coefficients of such functions which implies excellent approximation properties for the latter. Due to the finite numerical precision, rapidly decreasing functions can be seen as essentially periodic if the computational domain is chosen large enough such that the function and its first derivatives vanish at the boundaries of the domain. For such functions the Fourier coefficients decrease to machine precision (10 −16 in our case) for a moderate number N of Fourier modes. The two-dimensional Fourier transform is defined aŝ For gKP there is the additional advantage that the antiderivative is defined most conveniently in Fourier space. The Fourier multiplier −i/k x is regularized for k x → 0 as discussed in [15] by adding a small imaginary constant to k x . The gKP equation takes in Fourier space the form where L = ik 3 x − iλk 2 y /k x , and where F(u) = −ik x u n+1 /(n + 1). To avoid numerical errors due to the regularization of the term −i/k x , we rewrite (2.2) as an equation forŵ, whereû = ik xŵ , i.e., for the integral of u with respect to x. This is possible since we use an explicit scheme as explained below, and since the initial condition u 0 is chosen to be an x-derivative of a rapidly decreasing function. In addition the term F in (2.2) is an x-derivative itself, whereas L appears in the used scheme only in the exponential integrator. Thus we use (2.3) forŵ and F/(ik x ) in the actual computation which ensures that u is for all time steps an x-derivative within numerical precision.
An important practical advantage of Fourier transforms is the fact that the derivatives are diagonal in Fourier space. This allows for an efficient time integration since for equations of the form (2.2), there are many high-order time integrators, see e.g. [3,11,9,12,14] and references therein, especially for diagonal L as in the Fourier case. For the numerical integration of (2.2), the Fourier transform will be approximated via a truncated Fourier series, a discrete Fourier transform which will be computed with a fast Fourier transform. Thus the equations (1.1) will be approximated via a system of ordinary differential equations (ODEs) of finite dimension for the Fourier coefficients. The latter system is in the present case stiff, where the term stiffness is used to indicate that there are vastly different timescales in the studied problem. This makes the use of standard explicit methods inefficient for stability reasons.
It was shown in [14] that ETD schemes perform best for KP equations among the studied stiff integrators, and that the performance is similar for different ETD schemes. For these methods one uses a constant time step h = t m+1 − t m and integrates (2.2) with an exponential factor to obtainû The different ETD schemes differ in the approximation of the integral in (2.3). We use here the method by Cox and Matthews [3] of classical order four. An important aspect in the implementation of ETD schemes is the accurate and efficient computation of the so-called φ-functions, which appear in all ETD schemes. To avoid cancellation errors, we use here contour integrals in the complex plane as in [11] in the enhanced version [22] as discussed in [12].
Accuracy of the numerical solution is controlled as in [12,14] via the numerically computed mass (1.2) or energy (1.3) which will depend on time due to unavoidable numerical errors. We use the quantity (and similar for M ) as an indicator of the numerical precision. It was shown in [12,14] that the numerical accuracy of this quantity overestimates the L ∞ norm of the difference between numerical and exact solution by two to three orders of magnitude. A precondition for the usability of this quantity is a large enough number of Fourier modes. In practice the numerical error cannot be smaller than the smallest modulus of the Fourier coefficients. Note that the energy is more sensitive to a loss of regularity since it contains a derivative. Thus we generally consider this quantity unless otherwise noted. But the mass conservation is typically of the same order as the energy conservation in the numerical experiments.
We typically choose the number N of Fourier modes high enough that they decrease to machine precision for the initial data. The number N thus depends on the size of the computational domain, x ∈ [−π, π]L x , y ∈ [−π, π]L y where the real constants L x , L y are chosen large enough to ensure 'periodicity' of the initial data in the sense discussed above. For gKP there is the problem of the algebraic fall off of the solution for t > 0 even if the initial data are rapidly decreasing. Thus we are forced to use a larger domain than would be necessary for corresponding initial data for gKdV where the solution stays rapidly decreasing if the initial data are. In practice we reach a resolution in Fourier space of better than 10 −10 for times much smaller than the time of a possible blow-up, see also [14].
The occurrence of a blow-up leads to an increase of the Fourier modes for the high wave numbers which eventually breaks the code. Since we use an explicit method, we can afford sufficiently small time steps to have satisfactory resolution in time up to t ∼ t * . Typically we run out of resolution in Fourier space first. A blow-up is also identified via diverging norms of the solution (see below).

Dynamic rescaling.
In [13] we have used a dynamic rescaling of the gKdV equation to analyze blow-up in more detail with an adaptive approach. This method can be also generalized to gKP equations. The basic idea is to use the scaling invariance of the gKP equation discussed in the previous section, but now with a time dependent scaling factor. As for gKdV we consider the coordinate change with a factor L(t) This leads for (1.1) to where the index τ denotes the derivative with respect to τ . Thus the space and time scales are changed adaptively around blow-up which is reached here for τ → ∞. The asymmetry in x and y with respect to the rescaling with L also explains the stronger divergence of the y-derivative of u at blow-up than of the x-derivative as observed in [16].
As for gKdV, equation (2.6) is also important for theoretical purposes in describing asymptotically blow-up. In the limit τ → ∞, the functions U , v ξ , v η and a are expected to become independent of τ which is denoted by a superscript ∞. Thus (2.6) reduces in this limit to In contrast to gKdV, one does not get an ODE in this case, and there is no reason to assume that this partial differential equation reduces to an ODE in generic cases. There are in principle two different scenarios important in this context, an algebraic or an exponential decay of the scaling factor L(τ ). In the algebraic case, we have L(τ ) = C 1 τ γ1 with constants C 1 , γ 1 < −1/3 and thus a ∞ = 0 as well as Then equation (2.8) reduces for v ∞ η = 0 to the equation for travelling wave solutions of the gKP equations in a commoving frame which has the unique nontrivial localized solution Q (1.6). Note that the latter condition is automatically satisfied for initial data with a symmetry with respect to y → −y as in the examples we consider here. Since the gKP equation is invariant under this transformation, we have v η = 0 in the studied examples.
For exponential decay we have L(τ ) = C 2 e a ∞ τ with C 2 = const and a ∞ < 0. Relation (2.5) implies in this case For the numerical implementation, the scaling factor L and the speeds v ξ , v η have to be chosen in a convenient way. A possible choice for the latter is to fix the single (by assumption) global minimum of U at ξ = η = 0 which implies U 0 ξ = U 0 η = 0, where the superscript 0 denotes that the function is taken for ξ = η = 0. These conditions lead to In the numerical examples we study in this paper, it turns out that x m does not change much, whereas y m is constant for the reasons explained above. This is in contrast to gKdV in the L 2 critical case where it was proven in [20] that x m → ∞ for τ → ∞. Since the computation of the derivatives in (2.11) at each stage is expensive in two dimensions, we do not fix x m and y m in the shown examples.
The coordinate transformation (2.6) implies for the L 2 norm of u Thus the L 2 critical case is n = 4/3 as already mentioned. The L 2 norm of u x scales as which implies that it is invariant under the transformation (2.5) for n = 4. Since the blow-up theorems in [21] are established for the L 2 norm of u y , we consider here, (2.14) Fixing ||U η || to be constant, we get This will be chosen for the numerical implementation. The quantity L(τ ) and the physical time t are computed as in [13] via the trapezoidal rule. The accuracy of the numerical solution is controlled via the computed L 2 norm of U via (2.12) and the energy Note that the energy is invariant under the rescaling (2.5) as the L 2 norm of u x for n = 4. In this sense the case n = 4 is energy critical.
The spatial dependence in equation (2.6) as well as the time dependence will be treated as for the direct integration of gKP outlined in the previous subsection. The numerically problematic terms in (2.6) for the Fourier approach are ξU ξ and ηU η . In [13] the dispersive oscillations with slow decrease towards infinity caused for gKdV numerical errors at the boundaries via a pollution of the Fourier coefficients at the high wave numbers. The problem could be solved in this case by using a very high resolution in time to essentially propagate the solution with machine precision. But in two spatial dimensions, such a high time resolution is computationally expensive. For gKP there is the additional problem of the algebraic tails due to the antiderivative. The effects of these tails on the Fourier coefficients are much worse than those of the oscillations of small amplitude. As will be shown below, they make it impossible to compute long enough with the rescaled codes to clearly identify blow-up. Nonetheless the rescaling is important since we can use it via a postprocessing of the directly integrated solutions to identify the function L(t) from computed norms of the solution.

Numerical tests.
In addition we can compare the results of both codes which provides an useful test since the codes are independent. To this end we run the rescaled code for the examples with blow-up in the following sections as long as the numerically computed mass is conserved to better than 10%. Then we directly integrate gKP for the same initial data to the physical time corresponding to the end value of τ for the rescaled solution. The resulting solutions for the case n = 4/3 studied in Fig. 3.6 are shown on the x-axis in Fig. 2.1. We present both solutions in one figure on the x-axis where they differ the most (the solution to (2.6) is rescaled back to u). Both codes are run with N x = N y = 2 10 Fourier modes, and L x = L y = 5 for the direct integration (the values for the rescaled code are given in the following sections). It can be seen that the agreement is much better than indicated by the conservation of the numerically computed mass. The disagreement is mainly due to the imposed periodicity which poses a problem for the rescaled code since the dispersive radiation reenters the dynamically rescaled domain from the other side. In the same way we study the case n = 2 presented in Fig. 4.6. Again the main disagreement is due to the dispersive oscillations reentering the domain from the right. But the good agreement between both codes except for these oscillations provides the wanted test of the numerical approaches.
3. The L 2 critical case n = 4/3. In this section we will study solutions to the gKP equations with the L 2 critical nonlinearity n = 4/3. This nonlinearity is not very relevant for applications, but mathematically interesting since it corresponds to the gKdV case n = 4 which was intensively studied in [20] and references therein. A similar theoretical description could be possible in this case. In addition it is different from the supercritical cases n > 4/3, and thus will be treated in more detail. The third root of u is computed in standard way as u 1/3 = sign(u)|u| 1/3 to ensure that the real branch of the root is taken.
We first study for gKP I initial data with positive energy, namely u 0 = ∂ xx exp(−x 2 − y 2 ). The computation is carried out with L x = 20, L y = 4, N x = N y = 2 10 and N t = 1000 time steps for  (1.1) in blue on the left; on the right the same setting for n = 2 and the initial data U 0 = 6 ∂ ξξ exp(−(ξ 2 + η 2 )) at τ = 0.5. 0 < t ≤ 0.5. The relative computed energy is conserved to better than 10 −5 . As can be seen in Fig. 3.1, oscillations propagating to the left form, and the initial hump appears to be just radiated away. Due to the imposed periodicity, these oscillation reenter on the right. The algebraic tails to the right can be also clearly recognized. A consequence of these tails is that a slight Gibbs phenomenon appears which can be seen from the Fourier coefficients in the same figure. The solution is well resolved in y-direction, but the Fourier coefficients in k x no longer decrease to machine precision (but 10 −5 is more than sufficient for our purposes). There is no indication of blow-up, and this is even more obvious from the norms shown in Fig. 3

.2.
Both ||u|| ∞ and ||u y || 2 appear to be monotonically decreasing. This suggests again that the initial hump will be just radiated away towards infinity.
The situation is quite different for initial data with negative energy, u 0 = 12 ∂ xx exp(−x 2 − y 2 ). The calculation is performed with L x = L y = 5, N x = 2 10 , N y = 2 13 , and N t = 50000 time steps for 0 ≤ t < 0.078. As can be seen in Fig. 3.3, there are dispersive oscillations forming immediately and propagating to the left. Due to the imposed periodicity condition, these oscillations reenter at a given time on the right (note that only part of the computational domain is shown in the figure). But more importantly one can see that the initial minimum gets more and more peaked and finally appears to blow up in a point. The code is stopped once ∆ > 10 −3 .  The blow-up is even more obvious from the norms ||u|| ∞ and ||u y || 2 shown in Fig. 3.4 which clearly seem to diverge. The Fourier coefficients of the solution in Fig. 3.3 at the last recorded time are shown in Fig. 3.9. As expected from the rescaling (2.5), near a blow-up the y-derivative diverges more strongly than the x-derivative of the solution. This implies that the Fourier coefficients decrease much more slowly in k y than in k x . This is why a higher resolution in k y was chosen from the beginning. But this behavior will persist with even higher resolution in k y due the scaling of the solution close to blow-up. On the other hand the resolution in time appears to be sufficient since the results for half the number of time steps used here do not differ from the ones shown within the numerical limits imposed by the Fourier coefficients.  To understand the type of this blow-up better, we try to solve the rescaled equation (2.6) for these initial data. With L ξ = 11, L η = 10, N ξ = N η = 2 10 and N τ = 10 4 time steps for 0 < τ ≤ 0.1, we can solve (2.6) with a relative conservation of the computed mass of the order of 10 −1 . The solution at the final time can be seen in Fig. 3.6. It can be recognized that the dispersive oscillations and even more so the algebraic fall off of the solution lead to a Gibbs phenomenon at the boundary of the computational domain. This is reflected in the Fourier coefficients in Fig. 3.6, and by the fact that energy conservation is no longer given in this case. The problematic aspect is here the increase in the amplitude of the high wave numbers in k η . These lead eventually to a breakdown of the solution. Thus we could not study the solution on larger domains with more Fourier modes since even with a high resolution in τ , the code would crash.
The function a from (2.7) and the physical time in dependence of t can be seen in Fig. 3.7. Obviously we do not get close enough to the blow-up to decide whether a tends to zero or not. The physical time at the end of the computation also shows that we are not as close as necessary to t * to use this approach to discuss the type of blow-up. Therefore it is more promising to integrate gKP directly, and to infer the type of blow-up from the norms of the solution. If we assume that for n = 4/3 the asymptotic behavior (2.9) is given, we get with (2.5) as well as (2.14) (recall that γ 1 = −1 for gKdV in the case n = 4) We fit ln ||u|| ∞ for the solution shown in Fig. 3.3 to C 1 +c 1 ln(t * −t) and ln ||u y || 2 2 to C 2 +c 2 ln(t * −t) via the optimization algorithm [18] distributed with Matlab as fminsearch. Since the L 2 norm involves an integral over the whole computational domain, it is numerically more stable, and the fitting should be thus reliable. However, this is not the case if strong gradients appear in the solution as in the case of blow-up. It turns out that the L ∞ norm typically produces results in better accordance with the theory exactly since it only takes into account effects close to singularity. We present the results for both the L 2 norm of u y and the L ∞ norm of the solution to test the consistency of the result. The fitting is shown in Fig. 3.8 for the last 1000 computed points. We find for ||u y || 2 2 the values t * = 0.0763, c 1 = −3.1162 and C 1 = −1.4687, and for ||u|| ∞ we get t * = 0.0765, c 2 = −0.6751 and C 2 = 1.1063. The good agreement of the blow-up time t * shows the consistency of the approaches. The oscillations in the L ∞ norm suggest as already mentioned that the L 2 norm of u y should be more reliable. Note, however, that the fitted values for c 1 and c 2 are with (3.1) compatible with γ 1 = −1 (i.e., to −2 respectively −3/4), which corresponds to the gKdV case for n = 4 where L ∝ 1/τ . But the agreement is better for the L ∞ norm of u, presumable due to a loss of accuracy in the gradients at the boundaries of the computational domain which affect ||u y || 2 , but not ||u|| ∞ . Nonetheless the L 2 critical case seems to be for both gKdV and gKP characterized by an algebraic in τ vanishing of L.  Fig. 3.9. This shows that x m in fact changes with t, but much less so than in the L 2 critical case for gKdV where blow-up happens at infinity for almost solitonic initial data. Since we cannot get arbitrarily close to the actual blow-up for obvious reasons, it is difficult decide whether there will be blow-up at a finite x * . As already mentioned, y m = 0 for symmetry reasons. The plateaus in the plot Fig. 3.9 are due to the minimum just being evaluated on grid points. We fit ln x m as the norms of u above for the last 1000 time steps to α 1 ln(t * − t) + α 2 , where t * is one of the values determined above by fitting the norms of u. Doing this for the values obtained from the L ∞ norm of u, we get the figure shown in Fig. 3.9 and the values α 1 = −0.0974 and α 2 = 0.1569. Doing the same fitting for the t * obtained by fitting the L 2 norm of u y , we get α 1 = −0.0735 and α 2 = 0.3000. Thus the values do not coincide, but they are both negative which would imply a blow-up at infinity. As in the gKdV case [13], the asymptotic regime for the x m appears much later in the computations than the one for the norms of u, and as there, it is difficult to make reliable predictions about the value of x * . From the present computations, it cannot be ruled out that the blow-up in the L 2 critical case happens for infinite x * . But the small value of |α 1 | (much smaller than the critical exponents for the norms discussed above) makes a finite x * more probable. Recall that for gKdV in the L 2 critical case x m ∝ L, see [20].

The location of the global minimum is traced as a function of time in
Another interesting question is which initial data are simply radiated away and which blow-up. Since no stable solitons are known, it is not surprising that no stable structures appear in the computations. If we consider initial data as above of the form (1.5), we find that there is no blow-up for values of β ≤ 6, but that there is blow-up for β ≥ 7. Numerically it is difficult to decide at which value of β blow-up appears, since it is difficult to distinguish a strongly peaked minimum from an exploding norm of the solution. But the important observation is that the energy for initial data with β = 7 is positive. Thus negative energy does not appear to be the necessary condition for blow-up. But the mass of the initial data has to be large enough, and the energy has to be small enough for a blow-up to appear. Numerically we can only give indications where for certain classes of initial data the transition between global existence of the solution in time and finite time blow-up appears.
There are no theoretical results on blow-up for gKP II solutions. In [16], numerical results indicate that there is no blow-up in this context for n ≤ 2. We confirm this here by looking at solutions for n = 4/3 and u 0 = 6 ∂ xx exp(−(x 2 + y 2 )) with L x = 20, L y = 4, N x = 2 11 , N y = 2 10 and N t = 10 3 for t ≤ 2. We obtain a ∆ ∼ 10 −3 . As can be seen in Fig. 3.10, there is no indication of blow-up, the initial data seem to be simply radiated away. Note that the radiation propagates in this case to the right, whereas the algebraic tails stretch to −∞. Since we compute for a long time, the periodicity of the considered situation leads to a reentering of the oscillations from the left side. The norms for this solution in Fig. 3.11 indicate that the L ∞ norm decreases monotonically, whereas the L 2 norm of u y appears to reach a constant value, at least on the considered time scales. The latter seems to be due to Gibbs phenomena at the boundaries because of the algebraic decrease of the solution towards infinity. The decrease of the L ∞ norm of the solution shows in any case a completely different behavior from the gKP I blow-up scenarios. Note that the solutions behave qualitatively the same for twice the initial data in (3.10). Thus we did not find an indication for blow-up in gKP II solutions for n = 4/3. It seems that the defocusing character stabilizes the solution against blow-up, whereas the opposite is true for gKP I, which has a focusing effect. 4. The supercritical case n = 2. In this section we study the supercritical case n = 2 which appears to be typical for n > 4/3 for gKP I solutions. It is also relevant in applications since it appears in the modelling of sound waves in antiferromagnetic systems [6]. We will again consider the initial data (1.5) for various values of the parameter β.
First we study gKP I (1.1) with n = 2 and the initial data u 0 = ∂ xx exp(−(x 2 + y 2 )) for which the energy is positive. We use L x = 10, L y = 4, N x = N y = 2 10 and N t = 10 4 for t ≤ 0.1 and obtain a relative conservation of the computed energy of the order of 10 −8 . It can be seen in Fig. 4.1 that the initial data appear to be simply radiated away. The situation is completely different for initial data u 0 (x, y) = 6 ∂ xx exp −(x 2 + y 2 ) for which the energy is negative. The computation is carried out with L x = L y = 5, N x = 2 11 , N y = 2 13 and N t = 50000 time steps for t ∈ [0, 0.0265]. As can be seen in Fig. 4.3, the solution develops the same dispersive oscillations propagating to the left as in Fig. 4.1 (only part of the computational domain is shown), but this cannot stop the negative peak from becoming more and more compressed. This peak finally appears to blow up in a point. The code is stopped at the time t = 0.0258375 when ∆ > 10 −3 . Note that the final stage before blow-up is happening on much shorter time scales than in the previous section for n = 4/3.
It can be seen in Fig. 4.4 that we run out of resolution in Fourier space. As expected from the rescaling (2.5), there are strong gradients especially in y which cannot be addressed with even higher resolution.
Both the L 2 norm of u y and the L ∞ norm of u indicate a blow-up as is clear from Fig. 4.5.
In an attempt to identify the type of blow-up, we solve the rescaled gKP equation (2.6) for the initial data U 0 = 6 ∂ ξξ exp(−(ξ 2 + η 2 )). The computation is carried out with L ξ = 3, L η = 7, N ξ = N η = 2 10 and N τ = 10 4 for τ ≤ 0.5 resulting in a relative conservation of the computed mass of the order of 10 −3 (again energy conservation is much worse due to Gibbs phenomena at the boundary of the computational domain). The solution at the final time is shown in Fig. 4.6. It can be seen that the 'zooming in' effect of the rescaling close to blow-up is at work. But due to the dispersive oscillations with a slowly decaying amplitude and due to the algebraic fall off of the solution, this eventually leads to a Gibbs phenomenon. This is clearly reflected in the Fourier coefficients of the solution in the same figure. The instabilities due to the increase in the high wave numbers will eventually break the code. This problem does not disappear on a larger computational domain, the code tends to become unreliable even earlier due the above described problems.
In Fig. 4.7 we show the logarithmic derivative a in (2.7) of the scaling factor L. It decreases rapidly to strongly negative values which implies that most of the dynamical evolution happens for comparatively small τ . The quantity a appears to approach a negative constant which would indicate an exponential τ -dependence of L. The oscillations in a are due to the imposed periodicity forcing the dispersive radiation to reenter the computational domain on the right side. Still we do not get close enough to the blow-up to be able to decide which type of blow-up is realized here. This is also clear from the physical time we show in the same figure which is not close enough to the time in Fig. 4.3 where the direct integration breaks.   Due to the described instabilities of the code we did not succeed to run it with higher resolutions on larger domains even with a very high time resolution. Thus it appears again more promising to trace the norms of the solution obtained via a direct integration of the gKP I equation as in Fig. 4.3. The norms indicate as Fig. 4.7 an exponential τ dependence of L which implies with (2.5) and (2.14)  Solution to the rescaled gKP I equation (2.6) with n = 2 for the initial data U 0 = 6 ∂ ξξ exp(−(ξ 2 +η 2 )) at τ = 0.5 on the left, and the modulus of the corresponding Fourier coefficients on the right.
As before we fit ln ||u|| ∞ for the solution shown in Fig. 4.3 to ln C 1 + c 1 ln(t * − t) and ln ||u y || 2 2 to ln C 2 + c 2 ln(t * − t). Since strong gradients appear close to blow-up, the L ∞ norm of the solution appears here to be asymptotically a more reliable indicator. The fitting is performed for times close to the last recorded time since the asymptotic behavior can be only expected for t ∼ t * . To have enough points for a reliable fitting, the time interval is chosen to include the 800 last computed points up to t = 0.0257845. We get for ||u y || 2 2 the fitting parameters C 2 = −2.5795, c 2 = −0.9851 and t * = 0.0258. Thus c 2 is compatible with the theoretically expected −1. The quality of the fitting can be seen in Fig. 4.8. For ||u|| ∞ we get C 1 = 0.2878, c 1 = −0.4445 and t * = 0.0258. Despite the oscillations in ||u|| ∞ in Fig. 4.8, the value of c 1 is very close to the expected −1/3, and the value for t * coincides within numerical precision with the value found from ||u y || 2 which shows the consistency of the approach.
The location of the global minimum can be seen as a function of time in Fig. 4.9. The behavior is very similar to the case n = 4/3 in Fig. 3.9. Again it is difficult to decide whether there will be a blow-up at finite x * , whereas again y m = 0 for symmetry reasons. If we fit ln x m as the norms of u above for the last 800 time steps to α 1 ln(t * − t) + α 2 , where t * is one of the values determined above by fitting the norms of u, we get the figure shown in Fig. 4.9 and the values α 1 = −0.0488 and α 2 = 0.2281. Thus the negative value for α 1 could again imply a blow-up, but it is probable that we simply did not get close enough to the blow-up to be able to decide. The small value of |α 1 | is also consistent with a finite value of x * .  The above example clearly indicates a blow-up for solutions to the gKP I equation for n = 2 for initial data with sufficiently small energy. An interesting question is again to identify the condition on the initial data for blow up. For initial data of the form (1.5) we find that the solution stays regular for β ≤ 2. But there appears to be blow-up for initial data with β ≥ 3 for which the energy is positive. Again it is numerically difficult to get closer to the actual threshold.
Blow-up does not seem to appear though for the corresponding gKP II equation. If we consider the same initial data that led to blow-up for gKP I, the energy is again negative, but there is no indication of blow-up in this case as can be seen in Fig. 4.10. The computation is carried out with L x = 10, L y = 4, N x = N y = 2 10 and N t = 10 3 time steps for t < 0.1. The typical tails, for gKP II to the left, can be clearly seen as well as the dispersive oscillations propagating in the same direction. Because of the periodicity, they reenter on the right. But it appears that generic localized initial data will be just radiated away to infinity for gKP II for all n.
This is even more obvious from certain norms of the solution as shown in Fig. 4.11. Both the L 2 norm of u y and the L ∞ norm of u appear to decrease monotonically. Note that the same qualitative behavior is obtained for twice the initial data considered in Fig. 4.10. Thus we did not find an indication for blow-up in gKP II solutions for n = 2.  exclude blow-up for the latter, it definitely does not appear at comparable energies or masses as for gKP I, presumably due to the defocusing effect of gKP II. In this section we will study whether there is blow-up in gKP II solutions for n = 3 and n = 4. Since the behavior of gKP I solutions appears to be similar to the case n = 2, we concentrate here on gKP II solutions.
For initial data of small enough mass, the gKP II solutions appear again to be simply radiated away. But for initial data of the form (1.5), the situation changes for large enough β. In Fig. 5.1 we show the gKP II solution for β = 6 for several values of t. The computation is carried out with L x = 5, L y = 4 and N x = N y = 2 12 Fourier modes and N t = 20000 time steps for t < 0.0014. The code is stopped once the numerically computed relative energy is conserved to less than 10 −3 . It can be seen that this time the maximum of the solution appears to blow up.
The resolution in Fourier space at the last recorded time can be seen in Fig. 5.2. Again due to the behavior indicated by the rescaling (2.5), there are strong gradients especially in y which imply that the loss in resolution in Fourier space is predominantly in y-direction.  As for the gKP I blow-ups in the previous sections, both the L 2 norm of u y and the L ∞ norm of u indicate a blow-up as is clear from Fig. 5.3.
This already indicates a similar type of blow-up as for gKP I. This is confirmed by an asymptotic analysis of the norms in Fig. 5.3 close to the blow-up. As before we fit ln ||u|| ∞ for the solution shown in Fig. 5.1 to ln C 1 +c 1 ln(t * −t) and ln ||u y || 2 to ln C 2 +c 2 ln(t * −t). For the fitting of ||u|| ∞ we use the 500 last computed points and get C 1 = 1.8650, c 1 = −0.1721 and t * = 1.3335 * 10 −3 .  Despite the oscillations in ||u|| ∞ in Fig. 5.4, the value of c 1 is close to the expected −2/9. The fitting is less convincing for the L 2 norm of u y . For the last 200 time steps we get C 2 = −4.256, c 2 = −2.576 and t * = 1.365 * 10 −3 . The value of t * is clearly too large for what is to be expected by the breaking of the code. This indicates that the L 2 norm of u y is again affected by errors in the gradient close to the boundaries of the computational domain where we do not have enough resolution. But the results are compatible with what was obtained in the previous sections for gKP I, and the fitting for the L ∞ of u appears to be conclusive.
The same experiment will be carried for n = 4. For small enough β in (1.5), the initial data will again just be radiated away. This changes for β = 3 as can be seen in Fig. 5.5. The computation is carried out with L x = L y = 5 and N x = N y = 2 12 Fourier modes and N t = 20000 time steps for t < 0.0007. The code stops shortly after the numerically computed relative energy is conserved to less than 10 −3 , and we record only the data for which ∆ < 10 −3 . Since n is again even, the minimum blows up as before in these cases.
The resolution in Fourier space at the last recorded time can be seen in Fig. 5.6.
Both the L 2 norm of u y and the L ∞ norm of u clearly show a blow-up as can be seen in Fig. 5.7.
To identify the type of blow-up we again fit ln ||u|| ∞ for the solution shown in Fig. 5.5 to ln C 1 + c 1 ln(t * −t) and ln ||u y || 2 2 to ln C 2 +c 2 ln(t * −t). For the fitting of ||u|| ∞ we use the 500 last computed points and get C 1 = 0.91937, c 1 = −0.1623 and t * = 6.6953 * 10 −4 . The quality of the fitting can  be seen in Fig. 5.8. The value of c 1 is very close to the expected −1/6. The fitting is again less convincing for the L 2 norm of u y for the last 200 time steps. We get C 2 = 8.3069, c 2 = −0.858 and t * = 6.7126 * 10 −4 , which is, however, in accordance with expectations (c 1 = −2/3) and what was found for the L ∞ norm of u.
6. Conclusion. In this paper we have numerically solved the Cauchy problem for the gKP equations for smooth localized initial data with a single minimum and single maximum. The  results can be summarized in the following Conjecture: Consider the Cauchy problem for the gKP equations (1.1) with n ∈ Q, i.e., n = p/q with p, q ∈ N and p, q coprime, with initial data u 0 (x, y) with a single global minimum (for p even) or a single global maximum (for p odd) such that ∂ −1 x u ∈ S(R 2 ). Then • for n < 4/3, the solution is smooth for all t.
• for gKP II, the solution is smooth for all t for n ≤ 2.
• for gKP I with n = 4/3, initial data with sufficiently small energy and sufficiently large mass lead to blow-up at t * < ∞; asymptotically for t ∼ t * , the solution is given by a rescaled (via (2.5)) soliton (1.6) where the scaling factor L ∝ 1/τ for τ → ∞. This implies the blow-up is characterized by (6.1) • for gKP I with n > 4/3 and gKP II with n > 2, initial data with sufficiently small energy and sufficiently large mass lead to blow-up at t * < ∞; asymptotically for t ∼ t * , the solution is given by a localized solution to (2.8), which is conjectured to exist and to be unique, rescaled (via (2.5)) where the scaling factor L ∝ exp(κτ ) for τ → ∞ with κ a negative constant. This implies the blow-up is characterized by Obviously numerical experiments can just give indications on possible theorems, but the results presented in this paper are conclusive within the limitations discussed in the text.
An important open question is the condition on the initial data to lead to blow-up. It is unclear whether the relevant quantity is simply the energy, or whether this is related to the mass and energy of solitons in the cases where these exist (recall that there are no lumps for gKP I for n ≥ 4 and for gKP II in general). For gKdV it was proven in [20] in the L 2 critical case n = 4 that for initial data in the vicinity of the soliton, the criterion for blow-up is that the energy of the initial data is smaller than the soliton initial data, whereas the mass has to be greater than the soliton mass. Numerical experiments indicate, see for instance [13], that this is also the case for more general localized initial data. Thus it would be important to study solutions of (1.6), which are known to be unstable [4,5], and to see what type of perturbation if any leads to blow-up. To this end, these solitons, which are not known analytically, have to be constructed first. Note, however, that there is a whole family of soliton solutions parametrized by the speed c. For gKdV with n = 4 the mass of the solitons is independent of c, and the same holds for the energy which vanishes for all c.
Of similar interest as the solitons would be to solve the PDE (2.8) numerically which appears in the asymptotic description of blow-up and to check whether this gives the expected asymptotic description of the blow-up. A further interesting question is related to the asymptotic behavior of the location of the blow-up which presumably stays finite. But it is unclear how such a finite value is approached These questions will be the subject of further work.