Numerical study of shock formation in the dispersionless Kadomtsev-Petviashvili equation and dispersive regularizations††thanks: We thank B. Dubrovin, T. Grava, N. Tzvetkov and A. Weideman for helpful discussions and hints. This work has been supported by the project FroM-PDE funded by the European Research Council through the Advanced Investigator Grant Scheme, the ANR via the program ANR-09-BLAN-0117-01, the Austrian Science Foundation FWF, project SFB F41 (”VICOM”) and project I830-N13 (”LODIQUAS”). We are grateful for access to the HPC resources from GENCI-CINES/IDRIS (Grant 2013-106628) on which part of the computations in this paper has been done, and to the Vienna Scientific Cluster (VSC).
The formation of singularities in solutions to the dispersionless Kadomtsev-Petviashvili (dKP) equation is studied numerically for different classes of initial data. The asymptotic behavior of the Fourier coefficients is used to quantitatively identify the critical time and location and the type of the singularity. The approach is first tested in detail in dimensions for the known case of the Hopf equation, where it is shown that the break-up of the solution can be identified with prescribed accuracy. For dissipative regularizations of this shock formation as the Burgers’ equation and for dispersive regularizations as the Korteweg-de Vries equation, the Fourier coefficients indicate as expected global regularity of the solutions. The Kadomtsev-Petviashvili (KP) equation can be seen as a dispersive regularization of the dKP equation. The behavior of KP solutions for small dispersion parameter near a break-up of corresponding dKP solutions is studied. It is found that the difference between KP and dKP solutions for the same initial data at the critical point scales roughly as as for the Korteweg-de Vries equation.
Key words. Asymptotic Fourier analysis, Kadomtsev-Petviashvili equation, Burgers’ equation, Korteweg-de Vries equation, dispersive shocks
AMS subject classifications. Primary, 65M70; Secondary, 65L05, 65M20
Nonlinear evolution equations without dispersion and dissipation generically have solutions which show the wave breaking phenomenon, i.e., the formation of a shock, a gradient catastrophe in finite time. A standard example in this context is the Hopf equation
It is well known that the solution of (LABEL:Hopf) for an initial value problem can be obtained via the method of characteristics in the implicit form
If the initial data are such that is positive, the solution reaches a point of gradient catastrophe at where the derivative of the Hopf solution blows up, but where the solution stays finite.
Regularizations of this equation with small dissipation , the Burgers’ equation,
or small dispersion , the Korteweg-de Vries (KdV) equation,
will have solutions which stay regular at the shock of the Hopf solution for the same initial data, but show some critical behavior in the vicinity of the point . At the critical point, the difference between Hopf solution and the solution of the regularized equation shows a characteristic scaling in , for KdV . Dubrovin  conjectured a universal behavior of solutions to Hamiltonian regularizations of the Hopf equation (among which KdV is the most prominent example) in the vicinity of a shock. In [18, 19, 11, 20] strong numerical evidence for this conjecture was given, which was proven for KdV in  via Riemann-Hilbert techniques. For this difference scales as . An asymptotic description of dissipative regularizations was presented in .
Dubrovin’s conjecture is based on a double scaling limit where and simultaneously , in such a way that the limits
where , exist and are bounded. In [18, 19, 11, 20] it was shown that it is indeed possible to study these scalings in numerically if the critical point is known. The above formulae make it clear that the error in the values , must be smaller than the considered values of if the correct scalings are to be identified numerically. This implies that for a numerical study of critical scaling phenomena, a control of the error bounds for the critical point is crucial.
It is mathematically and physically interesting to study similar problems in higher dimensions where much less about shock formation and dissipative or dispersive regularizations is known. A breakdown of regularity always indicates a limit of the applicability of the studied model. It is at such points that effects as dissipation and dispersion, which have been neglected in the simplified model, become important. The main technical problem in higher dimensions is that the dispersionless system even of completely integrable equations as the Kadomtsev-Petviashvili (KP) equation , a -dimensional variant of KdV, is not integrable in the usual sense. This means that there are in general no standard solution generating techniques as hodograph methods or (linear) Riemann-Hilbert problems available in this context (but see  for a nonlinear Riemann-Hilbert approach). The KP equations read
where and where is a small scaling parameter. The case corresponds to the KP I model with a focusing effect, and the case corresponds to the KP II model with a defocusing effect.
The focus in this paper will be on the dispersionless variant of the KP equation, the dKP equation, also known as Khokhlov-Zabolotskaya equation , which follows from (LABEL:e1) for ,
It appears in many applications as a model for dissipation- and dispersionless, essentially one-dimensional waves with weak transverse effects, for instance in nonlinear acoustic and in gaz dynamics, see [46, 21]. It also plays a role in the context of general relativity and differential geometry [1, 12]. Alinhac and coworkers, see for instance  and references therein, showed that it appears as a universal model in the geometric study of singularity formation in nonlinear wave equations. Existence of solutions up to break-up was proven, and the singularity was identified for a generic setting as a one-dimensional cusp.
The dKP equation is not completely integrable in the sense that an infinite number of conserved quantities exists (in fact it only has three). But it is integrable according to the definition of  that an infinite number of hydrodynamical reductions exists, see also [26, 27, 47]. As was shown in , solutions can be constructed with Twistor methods in terms of Einstein-Weyl metrics. In [30, 31] solutions to the dKP equation were obtained via the solution of a nonlinear Riemann-Hilbert problem. This allowed the study of the long-time asymptotics and finite time break-up in [31, 28]. Another way to construct dKP solutions arises from the theory of Frobenius manifolds. In  the infinite dimensional Frobenius manifold corresponding to dKP was constructed. A common feature of the above solution generating techniques is that even the most explicit forms of the solutions require the solution of singular nonlinear integral equations and the inversion of implicit functions which makes it difficult to find explicit examples for the wave breaking. Since none of these approaches has been studied numerically so far, we solve here dKP directly up to the formation of a singularity.
As mentioned above, a precondition for the numerical study of scaling laws as for KdV is that the critical point and the critical solution can be determined with sufficient accuracy. It is the goal of this paper to provide the necessary tools for this in a rather general setting, and to apply them to the dKP equation for several classes of initial data. The applicability of this approach will be shown first for the -dimensional case where exact solutions provide tests. Note that a dissipative regularization does not give precise information on the critical point as can be inferred from the example of the Burgers’ equation with small dissipation we will also consider. The same will be shown to be true for a dispersive regularization as provided by KdV. Thus it can be concluded that a direct study of the dispersion- and dissipationless equation with the techniques explored in this paper is more promising in the context of identifying singularity formation.
The basic idea of the used numerical approach is to approximate the spatial dependence of the solutions by a discrete Fourier series. It is known that the asymptotic behavior (for large wave numbers ) of the Fourier coefficients is exponential for solutions analytic in the complex plane in the vicinity of the real axis. This dependence becomes algebraic at the break-up of the solution. Thus the vanishing of the exponential decrease in the Fourier coefficients would indicate the appearance of the singularity. In practice there are several problems in this context: firstly it is difficult numerically to integrate the equation up to the critical point; moreover even if this can be done with sufficient accuracy, the Fourier coefficients might be polluted especially at the high wave numbers where the asymptotic behavior is to be read off; and last but not least, the asymptotic relations for the Fourier coefficients are established for a continuous Fourier transform, whereas numerically only a discrete Fourier transform is considered as an approximation.
Asymptotic Fourier analysis was first applied numerically by Sulem, Sulem and Frisch , at the time of course with much lower resolution than is accessible today. Therefore the study in  was mainly qualitative, but gave a proof of concept. In  they studied singular solutions to the two-dimensional cubic NLS equation. An application of the method to the 2-d Euler equations can be found in [17, 32]. It has also been applied to the study of complex singularities of the 3-d Euler equations in , in thin jets with surface tension , the complex Burgers’ equation  and the Camassa-Holm equation  to name just a few examples. The main goal in these works is to decide whether the studied solutions develop a singularity. However, the purpose of the present paper is to identify time and location of a singularity, which is known to appear in the solution, with prescribed accuracy with these methods.
The tracking of singularities can be of course also achieved with other techniques, for instance Padé approximants in . The problem with the latter is, however, that this adds an additional dimension to the problem which is especially expensive if this technique is to be applied in higher dimensions. Full singularity tracking for a function in would imply the study of the function in , thus doubling the spatial dimensions. Therefore a quantitative numerical approach based on the asymptotic behavior of discrete Fourier series would be an economic way to study singularity formation in higher dimensions for functions analytic before a critical time . It is one goal of this paper to present a careful investigation of the possibilities and the limitations of this approach, and to show that the method can provide information on critical points with the needed accuracy to study scaling laws. The techniques are developped for the example of the Hopf equation and various regularizations thereof, and will then be applied to the dKP equation for various initial data. The thus identified dKP solution at the critical point is used to study the scaling in the small dispersion parameter of corresponding KP solutions. It is found that the scaling of the difference between KP and dKP solutions is compatible with . This suggests that the behaviour of the KP solutions in the singular direction is similar to KdV solutions near the break-up of Hopf solutions.
To establish asymptotic Fourier analysis as a quantitative tool to identify shock formation, we test the various numerical problems in the approach  separately and show they can be controlled. Then we will establish that this method in fact can be used to determine the critical point and the critical solution to the Hopf equation. The paper is organized as follows: In Sec. 2 we collect some known facts about the asymptotic behavior of Fourier coefficients of a real function analytic in the vicinity of the real axis. We numerically integrate the Hopf equation up to the critical point for given initial data and compare it with the exact solution. For the latter we show how well the Fourier coefficients can be fitted to the expected asymptotic decrease. We perform this fitting during the numerical solution of the Hopf equation to determine the critical time. The same approach is applied to solutions to the Burgers’s and KdV equation for the same initial data. In Sec. 3 we apply these techniques to the break-up in dKP solutions for two classes of initial data, localized in two dimensions or as the line solitons of KP, localized in one spatial direction and infinitely extended in the other. In both cases, the scaling of the difference between KP solutions in the small dispersion limit and dKP solutions at the critical point is studied. We add some concluding remarks in Sec. 4.
2 1+1-dimensional case
In this section we will address numerically -dimensional examples from the family of Hopf equations. We first review some well known facts from asymptotic Fourier analysis. Then we will numerically integrate the Hopf equation for a concrete example up to shock formation. We discuss how the asymptotic behaviour of the Fourier coefficients can be used to identify the critical time and solution of the Hopf equation. A similar analysis is then performed for the Burgers’ and the Korteweg-de Vries equation which can be seen as dissipative and respectively dispersive regularizations of the Hopf equation.
2.1 Asymptotic Fourier analysis
In this subsection we review some well known facts on the asymptotic behavior of the Fourier transform of a real function which is analytic in a strip around the real axis. The resulting formulae will be then used to numerically track singularities of the solutions in the complex plane via their Fourier coefficients without having to deal with two real dimensions as in . This allows to determine when a singularity hits the real axis, i.e., when the real solution becomes singular.
Spectral methods are a powerful tool in the numerical solution of differential equations because of their excellent approximation properties for smooth functions, see for instance . It is well known that the error in the approximation of such a function with a spectral approximant of order , e.g. a polynomial of degree , decreases faster than any power of . This provides a huge advantage over methods with just an algebraic decrease of the error and is the basis of the efficiency of spectral methods, see also . The most used such methods are the expansion of the studied function in terms of discrete Fourier series. For given analyticity properties of the studied function, precise mathematical statements on the asymptotic behavior of the Fourier coefficients exist. As we will detail below, these imply that the location of singularities in the complex plane can be obtained from a given Fourier series computed on the real axis. It should be possible in this way to numerically determine the type of the singularity, to trace it in the complex plane, and to establish when the singularity hits the real axis as was done for the first time in .
The Fourier transform of is defined as
A singularity in the complex plane of the form , , with in the lower half plane () implies with a steepest descent argument for the following asymptotic behavior of the Fourier coefficients (for a detailed derivation see e.g. ),
Consequently for a single such singularity with positive , the modulus of the Fourier coefficients decreases exponentially for large . For , i.e., a singularity on the real axis, the modulus of the Fourier coefficients has an algebraic dependence on . If there are several singularities of this form at , , there will be oscillations in the modulus of the Fourier coefficients for moderately large .
To numerically compute a Fourier transform, it has to be approximated by a discrete Fourier series which can be done efficiently via a fast Fourier transform (FFT), see e.g. . The discrete Fourier transform of the vector with components , where , (i.e., the Fourier transform on the interval where is a positive real number) will be always denoted by in the following. There is no obvious analogue of relation (LABEL:fourierasym) for a discrete Fourier series, but it can be seen as an approximation of the latter, which is also the basis of the numerical approach in the solution of the PDE. It is possible to establish bounds for the discrete series, see for instance . We will in this paper always fit the discrete Fourier series for large to the asymptotic formula (LABEL:fourierasym).
Note that the minimal distance in Fourier space is with being the number of Fourier modes and the length of the computational domain in physical space. This defines the smallest distance which can be resolved in Fourier space. All values of below this threshold cannot be distinguished numerically from 0.
2.2 Numerical integration of the Hopf equation
In this subsection we will study the integration of the Hopf equation for a concrete example up to the appearance of a gradient catastrophe. Since the method of characteristics gives an exact solution in this case in implicit form which can be obtained numerically in principle with machine precision, we can use this as a test case for the numerical integration of a nonlinear dispersionless equation up to the first critical point in the solution. We discuss how to obtain asymptotically reliable Fourier coefficients, and with which precision this can be achieved.
In the following we will consider the example for which the gradient catastrophe appears at the critical time and at the critical point . These initial data are motivated by the well known KdV soliton. They have the advantage that they belong to the class of rapidly decreasing functions which, for sufficiently large , can be analytically continued as periodic functions within numerical precision (the quantity and thus the length of the interval is chosen large enough that is smaller than machine precision ( in our case). This ensures that there is no Gibbs phenomenon and that the analytic behavior of the Fourier coefficients is not affected by effects at the boundary of the computational domain.
Numerically we solve equation (LABEL:Hopfsol) iteratively,
with the initial guess . If some relaxation is used, the iteration converges and is stopped once the residual of is smaller than in the norm. The iteration is done in a vectorized and thus efficient way, i.e., for all points at the same time. For the studied example the solution to the Hopf equation can be seen in Fig. LABEL:hopfexact at the initial and at the critical time.
To directly integrate the Hopf equation numerically, we use a standard Fourier spectral method in : for , we take , and the modulus of the Fourier coefficients computed via a discrete Fourier transformation decreases to machine precision for more than Fourier modes. Treating the spatial dependence in (LABEL:Hopf) with a discrete Fourier transform, we obtain for the Hopf equation a system of ordinary differential equations (ODEs) in for the Fourier coefficients. Since the latter system is not stiff in contrast to the KdV equation treated in a later subsection, we can use standard time integrators here. We choose the well-known explicit fourth order Runge-Kutta scheme for convenience.
The accuracy of the numerical solution is tested via comparison with the exact solution. In addition we trace a conserved quantity of the Hopf equation, the energy,
which, when computed for the numerical solution, will depend on time due to unavoidable numerical errors. As shown for instance in [23, 24], the conservation of energy can be used as an indicator of numerical accuracy. We test this for the Hopf equation below by tracing the quantity , where is the numerically computed energy.
Since the solution of the Hopf equation is very different for times and for , we will first numerically solve the Hopf equation up to . For Fourier modes, the solution at time is fully resolved as can be seen from Fig. LABEL:hopffourier_18, where the modulus of the Fourier coefficients is shown. The modulus goes down to machine precision. For a time step the norm of the difference between numerical solution and exact solution is , the quantity . For , the respective numbers are and , for one finds and . Thus overestimates the numerical precision by roughly two to three orders of magnitude (once it is of the order of machine precision, it is obviously dependent on rounding errors and can even increase with smaller resolution as above). It can be in fact used as an indicator of numerical accuracy in cases for which no exact solution is known, provided that there is sufficient spatial resolution. To sum up, the solution can be obtained up to machine precision for times .
For larger , time steps of different size might be convenient. Thus we run the code with a certain time step up to time and use the result as initial data for a second time evolution up to the critical time (since we have an exact solution for the studied example, it could be also used as exact initial data at that time; but since the goal is to develop and test methods for a general case for which no exact solution is known, this is not done here). The important point is that the solution can be obtained with machine precision at any time numerically.
For the studied example, one gets for with an norm of the difference between numerical and exact solution of the order of 0.032. This error (as well as the norm of the difference between numerical and exact solution) does not decrease for smaller time steps. This shows that the error, which is always biggest at the critical point as can be seen in Fig. LABEL:hopfc14, is due to a lack of resolution in and not in time. This lack of resolution is clear from the Fourier coefficients which can be also seen in Fig. LABEL:hopfc14. Visibly the Fourier coefficients for high wave numbers do not show the expected decrease with due to numerical errors at the critical time. This will be important for the asymptotic analysis in the following section.
Note that the quantity for the computed energy decreases to machine precision if the time step decreases to . This shows that this quantity is only a useful indicator of the numerical accuracy for sufficient spatial solution, i.e., if the Fourier coefficients decrease to the wanted accuracy. Lower values of than the modulus of the Fourier coefficient of the highest wave numvers are meaningless. To obtain higher accuracies, it appears necessary to allow for higher spatial resolution. For Fourier modes and , we get , for and we find . The time steps have been chosen in a way that the error does not decrease further if is decreased.
2.3 Asymptotic fitting of the Fourier coefficients
In this subsection we will study the asymptotic behavior of the Fourier coefficients for the solution to the Hopf equation discussed in the previous subsection. We show how to apply the asymptotic formula (LABEL:fourierasym) to the exact solution before and at the critical time. This will be done here first for the exact solution to avoid the problems with errors in the numerically determined Fourier coefficients at the critical time for high wave numbers.
As in the previous subsection, we discuss the Hopf equation (LABEL:Hopf) for the initial data . For a time , the solution and the Fourier coefficients computed via FFT with Fourier modes on the interval with can be seen in Fig. LABEL:hopfex18. The latter show the expected exponential decrease up to a level of roughly where the error saturates because of the finite numerical precision.
In this case, only one singularity is expected in the lower half plane. To test the asymptotic formula (LABEL:fourierasym) in practice, we first do a least square fitting for , which should be according to (LABEL:fourierasym) of the form
The fitting is done for a given range of wave numbers (we only consider positive ). For obvious reasons, has to be limited to values for which is larger than the rounding errors. We choose this threshold to be . In principle formula (LABEL:fourierasym) holds only for . However if we do the fitting for all positive , we get , and as well as , where we define the quantity as the norm of the difference between the solution and the the asymptotic formula (LABEL:abd),
As expected the difference is maximal near . For we get , and and for the quantity . The difference between and the fitted formula is shown in Fig. LABEL:hopfex18fit. The main difference appears for small . From the figure it can be concluded that a minimal value of might be optimal. On the other hand it can be also seen that rounding errors lead to some fuzzy behavior of the Fourier coefficients at high wave numbers. This suggests to consider only values . With these choices, we obtain , and . The difference between and the fitted formula can be seen in Fig. LABEL:hopfex18fit, it is of the order . Note that the value of is not much affected by the different choices of these thresholds. This just reflects the fact that the fitting to an exponential is much more stable than the fitting to an algebraic function.
The norm (LABEL:Delta) of the difference between and the fitted formula (LABEL:fourierasym) can be used to indicate the quality of the fitting. An additional test is provided from the relation between and . Since we compute the Fourier coefficients via a discrete Fourier transform whereas formula (LABEL:fourierasym) holds for a standard continuous Fourier transform, an additional term appears in the formula, . We get for the right hand side of this relation , a value close to the computed . However it is clear that the quantity is the more convenient indicator of the quality of the fitting. Thus it will be used to this end in the following.
The above fitting has only been carried out for the absolute value of the Fourier coefficients. To determine , the real part of the location of the singularity, we fit the imaginary part of the logarithm of ,
Since the logarithm is branched in Matlab at the negative real axis with jumps of , the computed will in general have many jumps. Thus one has first to construct a continuous function from the computed . The analytic continuation will be done in the following way: starting from the first value (largest ), we check for all other values of whether . If this is the case, we put . The result of this procedure will be a continuous function which will be fitted with a least square approach to a linear function. For the above example, we get and . The norm of the difference of the function and the fitted straight line,
can be seen in Fig. LABEL:hopfex18fita. It is of the same order of magnitude as the quantity which shows the self consistency of the approach.
The above situation is typical for a single singularity in the complex plane away from the real axis. This picture changes if , i.e., when the singularity hits the real axis which happens for the considered example at the critical time . In this case the Fourier coefficients as expected no longer show an exponential decrease as can be seen in Fig. LABEL:hopfexcfourier. Note the difference to the numerically computed Fourier coefficients at high wave numbers in Fig. LABEL:hopfc14.
The asymptotic fitting of as in Fig. LABEL:hopfex18fit for yields , and . The difference between and the fitting curve can be seen in Fig. LABEL:hopfexcfit. Since there is no exponential decay in this case, the fitting is much more sensitive to the considered cutoffs, for low due to the condition for the asymptotic formula to be valid, for large due to rounding errors. This is also obvious from Fig. LABEL:hopfexcfit where the errors are maximal at both ends of the shown spectrum. For we get , and with a fitting error as shown in Fig. LABEL:hopfexcfit. In both cases is of the order of . The quantity should be equal to and is as before much more sensitive to a choice of the limits and for the fitting.
By fitting the quantity for , we get and with a fitting error of the same order. The former is very close to the exact value . If we repeat the above fitting for Fourier modes for ( as before), we obtain , and with a and and . Thus a higher number of Fourier modes makes the fitting more stable and reliable.
To sum up, we have shown in this subsection that the quality of the fitting depends on the choice of the minimal and maximal values for the wave numbers. Generally a choice of or higher appears to be recommended. However this has the disadvantage that many Fourier modes are excluded, especially at times where only few modes contribute. Rounding errors enforce a choice of of the order of . By studying the difference between the considered quantity and the fitting curve, these values can be optimized as outlined in the next subsection. A proper choice of these thresholds is more important at the critical time when the singularity hits the real axis since there is no exponential decay in this case.
2.4 Asymptotic behaviour of the Fourier coefficients for a numerical solution to the Hopf equation
In the previous subsection we have studied the asymptotic fitting of the Fourier coefficients on a finite interval of the wave numbers and have fixed the limits of this interval essentially by hand. Here we will do this in a more systematic way by choosing a suitable lower threshold (which depends on the studied problem), and by varying the upper limit to reach a prescribed difference between the Fourier coefficients and the fitted curve with a maximal number of coefficients. Since this difference would obviously be minimal if upper and lower limit coincide, we impose the additional condition that at least half of the Fourier coefficients with will be included in the fitting to include most of the available information on the solution. With this approach and the preliminary studies of the previous sections, we are then able to perform the asymptotic fitting during the actual numerical computation of the Hopf solution. We show that the vanishing of the computed can be in fact used to determine the critical time and the type of the singularity via the exponent . The goal is to obtain an error in the fitting of the same order as the numerical error.
The fitting method used in this paper is obviously not the only possible one. For instance in  and related publications the sliding fit approach is applied where only few values of are taken into account in the procedure (typically between and for a resolution of 1024 modes, see for example ). The quality of the fitting in this approach is roughly determined by the independence of the computed fitting parameters on the minimal value of the used in the fitting, and on the discretization size . Note that the sliding procedure confirms the results found here, for and up to points considered. Nevertheless, this approach is computationally expensive which implies that the time of the singularity formation has to be previously roughly determined in another way to be able to perform the study close to it. Since it is the aim of this paper to identify the break-up numerically via the asymptotic behavior of the Fourier coefficients, even in cases where the critical time is not known, we do not use the sliding fit here.
As in the previous section, the asymptotic fitting of is done on the interval , where we put
where denotes a prescribed value for in (LABEL:Delta) on the fit interval. The choice allows the fitting even for short times where there are only few Fourier coefficients with . It satisfies the condition whilst maximizing the number of the numerically computed Fourier coefficients included in the fitting procedure. Indeed, for given , larger values of would lead to a too small number of coefficients used for the fitting.
For , we observe the influence of on the time of vanishing of and the corresponding value of in Fig. LABEL:HRKdap and in Table LABEL:tab1. Due to the large computational domain considered (, ), a sufficient resolution has to be used to get a sufficiently small value of (, for and ).
If decreases to , one finds that . But for smaller values of , for example for , the values for both and become worse approximations of the exact ones. In fact, to get such small fitting errors, we have to consider too few Fourier coefficients which explains the lack of precision for these. Thus we add the requirement that at least half of the total number of the Fourier coefficients available have to be used for the fitting. For , this implies that the ‘minimal’ achievable fitting error whilst using enough Fourier coefficients is , which corresponds roughly to the precision with which we compute the numerical solution, see Sec. LABEL:numhopf.
Thus for different values of , and the related precision with which we compute the numerical solution, we observe in Fig. LABEL:HRKp01 the time evolution of and , and give the values of and at two different times in Table LABEL:tab2.
As increases, and without further restrictions on . The time dependence of
the fitting parameters is very similar for different , which indicates that the fitting is reliable.
We conclude from this example that the method can be used in practice to determine the critical time and the type of the singularity via the value of if sufficient spatial resolution is used. The choice of the minimal wave number for the fitting, and the achievable minimal fitting error depend on the studied problem.
2.5 Burgers’ equation
The Burgers’ equation (LABEL:burgers) with small dissipation , can be viewed as a purely dissipative regularization of the Hopf equation. An asymptotic description at the critical point was given in . Since the system of ODEs resulting from the approximation of the solution via a discrete Fourier series in the spatial coordinate is now stiff, we use a stiff integrator, a fourth order exponential time differencing scheme by Cox and Matthews  as implemented in [39, 23]. This method is very efficient for both Burgers’ and KdV equation, see e.g. . In both cases we consider as before the initial data for .
For the Burgers’ equation the shock of the Hopf solution for will be regularized as can be seen for in Fig. LABEL:burgers01 at time . There is a zone with a strong, but finite gradient of order .
The Fourier coefficients for this solution are also shown in Fig. LABEL:burgers01. It can be seen that they decrease (essentially exponentially) to machine precision. Doing a least square fitting for the modulus of the Fourier coefficients with for (this is the result of the procedure outlined in the previous section with a fitting error of the order ), we find , and . This shows that the parameter indicating the distance between the nearest singularity to the real axis stays finite even for as expected. The difference between the modulus of the Fourier coefficients and the fitted asymptotic formula is of the order of . Fitting the quantity in (LABEL:phi), we get and . This implies that the real part of the location of the singularity is shifted towards larger values than for the critical value of the Hopf solution which is not surprising since the Hopf shock also propagates with . The difference between and the fitted asymptotic formula is of the order of , too.
For times the solution in the vicinity of the Hopf shock becomes steeper as shown in Fig. LABEL:burgers04, which means that more Fourier modes have to be used in order that the modulus of the coefficients decreases to machine precision ( instead of before).
A least square fitting as before for the modulus of the Fourier coefficients for gives , and . The difference between the modulus of the Fourier coefficients and the fitted asymptotic formula is of the order of . Fitting the quantity in (LABEL:phi), we get and with a similar quality of the fitting.
If we trace the quantities and during the computation as a function of time for (for there are not enough Fourier coefficients above the level of rounding errors for small times), we get the expected behavior as can be seen in Fig. LABEL:burgersalpha. The imaginary part of the singularity decreases strongly with before the critical time of the Hopf solution and comes close to the axis there. After the decrease becomes very slow, and the singularity stays at a finite distance from the real axis. This shows that the time cannot be really inferred from the Burgers’ equation with small dissipation in a reliable way. The transition in all characteristic quantities of the solution is not distinguished at this time, and even for the quantity continues to decrease. Therefore to identify and , it appears better to study the Fourier coefficients of the Hopf solution without dissipation to determine the critical point. The real part of the singularity can also be seen in Fig. LABEL:burgersalpha. It describes in some sense the motion of the dissipative shock.
2.6 Korteweg-de Vries equation
The KdV equation (LABEL:KdV) with small dispersion can be seen as a purely dispersive regularization of the Hopf equation. For the initial data studied there, there will be a zone of rapid modulated oscillations forming near the shock of the corresponding Hopf solution for greater than the critical time as can be seen for in Fig. LABEL:kdv1e4 at time .
The Fourier coefficients for this solution are also shown in Fig. LABEL:kdv1e4. It can be seen that they decrease (essentially exponentially) to machine precision. Doing a least square fitting as before for the modulus of the Fourier coefficients with for (no further restrictions), we find , and . This shows that the parameter indicating the distance between the nearest singularity to the real axis stays finite even for as expected. The difference between the modulus of the Fourier coefficients and the fitted asymptotic formula is shown in Fig. LABEL:kdv1e4fit. It can be seen that there are damped oscillations for small which die out for large and which indicate the presence of several singularities in the complex plane in the vicinity of the real axis. A fitting to more than one singularity is problematic given the limited number of Fourier coefficients. Our study is here of a more qualitative nature.
Fitting the quantity in (LABEL:phi), we get and . This implies that the real part of the location of the singularity is shifted towards larger values than for the critical value of the Hopf solution which is again not surprising since the Hopf shock propagates. The difference between and the fitted asymptotic formula is also shown in Fig. LABEL:kdv1e4fit. There are damped oscillations in this difference for small , too.
For times there are many oscillations in the oscillatory zone as shown in Fig. LABEL:kdv1e404. The higher number of oscillations also implies more oscillations in the Fourier coefficients as can be seen in Fig. LABEL:kdv1e404.
Consequently a least square fitting as before for the modulus of the Fourier coefficients for (to exclude the last big oscillation in the Fourier coefficients which is presumably a numerical artifact), we find , and . The difference between the modulus of the Fourier coefficients and the fitted asymptotic formula is shown in Fig. LABEL:kdv1e404fit. In contrast to the case shown in Fig. LABEL:kdv1e4fit where there are essentially just two big oscillations, there are many more here. For the modulus of the Fourier coefficients this implies with the asymptotic formula (LABEL:fourierasym) that one gets in the first case just standard damped oscillations of harmonic type. In the latter case, several singularities in the complex plane approach the real axis which leads to more complicated oscillations (several phases) as can be seen in Fig. LABEL:kdv1e404fit. Fitting the quantity in (LABEL:phi), we get and . The additional structure in the oscillations of the Fourier coefficients implies that a linear fitting of the phase is not efficient as can be seen in Fig. LABEL:kdv1e404fit.
If we trace the quantities and during the computation as a function of time for (again there are not enough non-trivial Fourier coefficients for at small times), we get the expected behavior as can be seen in Fig. LABEL:kdv1e4alpha. The used fitting just takes care of the singularity in the complex plane closest to the real axis. The imaginary part of the singularity decreases strongly with before the critical time of the Hopf solution and comes close to the axis there. After the decrease becomes very slow, and the singularity stays at a finite distance from the real axis. But again there is no clear indication of the precise values of the critical quantities. The real part of the singularity can also be seen in Fig. LABEL:kdv1e4alpha. It gives in some sense the almost constant speed of the dispersive shock.
3 Kadomtsev-Petviashvili equations
In this section we will study shock formation in the dKP equation with the methods discussed for the example of the Hopf equation in the previous section. This will be done for both the dKP I and dKP II equations for two classes of initial data. The first class is motivated by the line solitons of KP, i.e., solutions exponentially localized in one spatial dimension and infinitely extended in the other. The second class of initial data are localized in both spatial directions. Then we will study for several values of the small dispersion parameter how the difference between the dKP solution and the corresponding KP solutions scales at the critical point and for in dependence of .
3.1 Theoretical preliminaries
We fist collect some known analytic facts on KP and dKP equations we will need in the following. The KP equation (LABEL:e1) is not in the standard form for a Cauchy problem since is not a timelike coordinate if it is discussed as a standard second order PDE, see e.g. . To be able to treat initial value problems in , equation (LABEL:e1) can be rewritten in evolutionary form,
The antiderivative in (LABEL:kpevol) is to be understood as the Fourier multiplier with the singular symbol . In the numerical computation this multiplier is regularized as described in . Equations (LABEL:kpevol) and (LABEL:e1) are equivalent for certain classes of boundary conditions as periodic or rapidly decreasing at infinity, i.e., the ones studied in the present paper. The evolutionary form of dKP follows from (LABEL:kpevol) for .
The divergence structure of the KP equations (LABEL:e1) has the well known consequence that
even if this constraint is not verified at the initial time. As shown in [14, 33] the solution to a Cauchy problem not satisfying the constraint will not be smooth in time for . This leads to numerical problems which could have a negative impact on the determination of break-up singularities in dKP equations, see . Thus we consider only initial data subject to the constraint (LABEL:const). The first example with period in is
It is known (see e.g.  for references and examples) that solutions to KP equations for initial data in the Schwarzian space of rapidly decreasing functions generically do not stay in this space, but develop tails with an algebraic fall off to infinity. These tails will lead to a Gibbs phenomenon in a periodic setting. To reduce the latter, we will choose a larger computational domain than would be necessary if the solutions remained in the Schwarzian space as for instance in the case of KdV. We always give the Fourier coefficients to show that these algebraic tails do not affect the results in our examples.
The dKP equation has only three conserved quantities and thus does not belong to the family of completely integrable equations having an infinite number of conserved quantities. Therefore standard solution generating techniques as dressing transformations  and linear Riemann-Hilbert problems cannot be applied in this context. In  it was shown that solutions to the dKP equation can be found in terms of Einstein-Weyl geometries which can be constructed with Twistor methods, see for instance [35, 44]. However this approach is rather implicit if one wants to construct solutions for a given Cauchy problem and has not yet been used to this end. In  it was shown that the dKP equation allows for an infinite number of hydrodynamic reductions and is thus integrable in this sense. Alinhac and coworkers  found that the dKP equation appears as a universal model in the geometric analysis of blow-up in nonlinear wave equations. They showed that finite time wave breaking is generic for solutions to this equation, proved existence of regular solutions up to the break-up time and classified the singularity. Localized initial data become singular in a point and the type of singularity is a one-dimensional cusp as for the Hopf equation, the second spatial dimension remains regular. However these methods do not provide formulae to determine for given initial data.
similar to the solution of the Hopf equation via the method of characteristics. But is here an integral over one component of the solution of a nonlinear vector Riemann-Hilbert problem, see  for details. The latter is equivalent to a systems of two coupled nonlinear, singular integral equations. If the latter system of integral equations is solved for a given normalization condition, the dKP solution follows from the implicit relation (LABEL:uF). To get a numerical solution to dKP in this way is not straight forward and has not been achieved for general initial data so far. But the form of the solution allows the study of singularity formation which has been done in [28, 29]. It is shown that there will be generically a gradient catastrophe at in one spatial direction, whereas the solution remains smooth in the second spatial direction. Formulae are given for the solution for small , and , and the type of singularity is identified again as cubic in the singular direction. In particular, they show that at the breaking time , if the initial condition for the dispersionless KP equation is even in , then the first breaking point is on the -axis, i.e., , and that the -derivative of is zero at . Our numerical results for the case of localized initial data in Sect. 3.3. are in accordance with this prediction. But it is not clear how to obtain in the formalism  the values for the critical point for given initial data. Therefore we solve the dKP equation in this paper directly.
The above mentioned theoretical results in [2, 28] and references therein indicate that break-up in dKP solutions is one-dimensional. We choose initial data symmetric with respect to the -coordinate to ensure a blow-up of the -gradient only. However this is without loss of generality since the same approach as applied below can be used for general initial data. In this case one first has to determine numerically the direction of the strongest gradient which will indicate where break-up occurs. Since the rest is as for the symmetric initial data, we concentrate on this conceptionally simpler case.
3.2 Numerical approach
We use here again a Fourier spectral method for the spatial coordinates and as for the Burgers’ and KdV equations in the previous section a fourth order exponential time differencing scheme derived by Cox and Matthews in  for the time integration. We consider periodic (up to numerical precision) solutions in and , i.e. solutions on . The computations are carried out with points for .
As discussed in the previous section, for a reliable identification of the shock formation via the asymptotic behavior of the Fourier coefficients, sufficient spatial resolution is crucial. In practice we needed to modes for the Hopf equation in dimensions. To allow such high resolution simulations, we had to parellelize the codes also for the asymptotic analysis of the Fourier coefficients. A prerequisite for parallel numerical algorithms is that sufficient independent computations can be identified for each processor, that require only small amounts of data to be communicated between independent computations. To this end, we perform a data decomposition, which makes it possible to do basic operations on each object in the data domain (vector, matrix…) to be executed safely in parallel by the available processors.
Our domain decomposition is implemented by developing a code describing the local computations and local data structures for a single process. Global arrays are divided in the following way: denoting by , the respective discretizations of and in the corresponding computational domain, is then represented by an matrix. The latter is distributed among processors such that each processor , ( denoting the number of processors we can access) will receive elements of corresponding to the elements
in the global array, and then each parallel task works on a portion of the data.
While processors execute an operation, they may need values from other processors. The above domain decomposition has been chosen such that the distribution of operations is balanced and that the communication is minimized. The access to remote elements has been implemented via explicit communications, using sub-routines of the MPI library. The computation of the 2d-FFT is performed by using a transposition approach which takes advantage of the existing sequential FFT routines (more precisely, we use serial FFTW routines ). The asymptotic fitting of the Fourier coefficients, here in one spatial direction, requires in addition 2 local communications, a core being reserved to do the required computations. More precisely, the processor which ‘takes care’ of sends these data to the reserved processor, which receives the data (thus 2 local communications per time iteration, a SEND, and a RECEIVE), and then it alone performs the least square fit required, which takes less time than the computation of the solution for one iteration, performed by the other processors. The main technical problem in the use of ETD schemes is the efficient and accurate numerical evaluation of the functions
i.e., functions of the form and higher order generalizations thereof, as explained in  and in . We use the approach given in  in our implementation. The computation of these functions takes only negligible time for the -dimensional equations studied here, especially since it has to be done only once during the time evolution. In addition each processor only computes the portion of these functions that the further local computations require, according to the domain decomposition specified above.
To control the numerical accuracy we trace as for the Hopf equation in the previous section and as in  the conservation of the norm of , the mass of the solution. It is exactly conserved for the KP equations, but due to unavoidable numerical errors the computed mass will numerically depend on time. Since mass conservation is not implemented in the code, it provides a valid check of the numerical accuracy for sufficient resolution in Fourier space. It was shown in  that the quantity
typically overestimates the norm of the difference between numerical and exact solution by two orders of magnitude.
3.3 Shock formation in dKP solutions for initial data localized in one spatial dimension
We will first study numerically the appearence of break-up singularities in solutions to the dKP equations for the initial data (LABEL:uini2) with the methods explained above. Whereas initial data localized in both spatial dimensions will have a gradient catastrophe in just one point, initial data localized in one dimension and infinitely extended in the other appear to blow up on a curve.
To determine a solution for the initial data (LABEL:uini2), we use . We observe that the solution develops a shock in the -direction at time , as can be seen in Fig. LABEL:u1uts2d.
To determine the time of the appearance of this shock, we apply the fitting procedure described in the previous section for the Fourier coefficients (again denoted by ) in -direction, i.e., , to formula (LABEL:abd). As before (see Section 2.4) we use at least half of the Fourier coefficients with values above the rounding error, and an appropriately chosen interval . In this way, for we find the minimal possible precision to be , for and . For this prescribed error, we find that , where vanishes, and . We observe at this time the formation of a shock in the -direction, which can be seen in Fig. LABEL:u1uts, where we show the behavior of the solution at different times, plotted on the -axis, and the corresponding Fourier coefficients, plotted on the -axis.
The time dependence of the fitting parameters is shown in Fig. LABEL:u1f1.
As expected, we observed a rapid decrease of , and the change of from a value to the value when approaching the critical time, determined by the vanishing of . In Fig. LABEL:u1fitp we show the time dependence of the fitting parameters close to , for different values of . By doing a fitting of , we find that at , .
In this case, we reach an error in the fitting of , which is the same as obtained for the Hopf equation. Note also that the approximation of is of the same precision as in the latter. However, such a precision cannot always be achieved as will be seen below.
The time evolution of the quantity in (LABEL:deltaE), used as an indicator of the numerical accuracy, and of the -norm of until the critical time are shown in Fig. LABEL:massampl.
The quantity indicating the numerical mass conservation increases close to , but stays below , which implies that the system is still well resolved up to the shock formation. Whereas the norm of stays finite, there are strong indications of a blow-up of close to the critical time, which indicates that the fitting in fact identified correctly the gradient catastrophe.
As one can see in Fig. LABEL:derivativesofu, where and are plotted at , the gradient catastrophe does not appear in only one spatial point for the infinitely extended initial data.
The situation is rather similar for the dKP II case, for the same initial data (LABEL:uini2). By using the same fitting bounds as before, ( and ), we find that the solution of the dKP II equation develops a shock at , which is exactly the previously identified value for the critical time for dKP I. The precision we could achieve here is also similar, , and . The solution at is shown in Fig. LABEL:udkp2exp.
A noticeable difference appears however in the profile of . One can observe in Fig. LABEL:uxuytcdkp2exp, that in contrast to the dKP I case, the strongest value of the -gradient here occur at the trailing part of the solution, whereas in the case of dKP I, it was occurring at roughly . This is related to the focusing and defocusing character of dKP I and dKP II respectively. In the first case, the most advanced point in propagation direction gets focused, whereas in the latter case, this happens for the trailing points.
3.4 Shock formation in dKP solutions for initial data localized in both spatial dimensions
In this subsection we will study shock formation in dKP solutions for initial data of the form (LABEL:uini1) which are localized in both spatial directions. It is known, see [2, 28], that localized data with a single maximum will lead to a break-up in a single point as expected.
In order to satisfy the constraint (LABEL:const), the data studied here do not have a single maximum. Therefore there will be actually two points of gradient catastrophe. The second gradient catastrophe appears for dKP I for negative for a slightly larger (for dKP II the role of the two critical points is interchanged). Since the code can only be reliably run until the first shock formation, we cannot obtain the second break-up with the method, and we will only discuss this time in the following.
For and , we observe that the solution develops a shock in the -direction at a time , as we can see in Fig. LABEL:u2uts2d.
This is in accordance with the experiments in  for the same initial data, where the break-up was studied qualitatively for a dissipative regularization of the dKP equation. There it was found that the gradient catastrophe for is reached roughly at . In Fig. LABEL:u2uts2d it can also be seen that the initially rapidly decreasing function develops tails with an algebraic fall-off to infinity. Due to the imposed periodicity condition, this leads to a Gibbs phenomenon at the boundary of the computational domain. To reduce the effect of this phenomenon on the Fourier coefficients, we have used a considerably larger computational domain than shown in Fig. LABEL:u2uts2d. We present in Fig. LABEL:u2uts the solution to the dKP I equation with initial data (LABEL:uini1), plotted on the -axis for several values of and the corresponding Fourier coefficients, plotted on the -axis. From the latter it can be inferred that the asymptotic behavior of the Fourier coefficients is dominated by the break-up singularity and not by the algebraic fall off of the solution.
In this case, we find that the minimal attainable fitting error (with the above described requirement to use at least half of the Fourier coefficients) is for and . For such bounds, vanishes at , where and . The considerably larger error in the fitting than for the initial data (LABEL:uini2) is related to the appearence of a second break-up singularity for the data (LABEL:uini1) as explained in Remark LABEL:remark. Nonetheless the results of the fitting are reliable as can be seen in the following. A fitting of yields . The value of is compatible with the expected value . The time dependence of the fitting parameters is shown in Fig. LABEL:u2f1 for .
The time dependence of the quantity (LABEL:deltaE) controlling the numerical error, and of until are shown in Fig. LABEL:massamplu2.
In this case, the gradient catastrophe appears in only one point , as we can see from Fig. LABEL:derivativesofu2, where we show the behavior of (left) and of (right) at . Note that the gradient in -direction is much smaller than in -directions which confirms the theoretical expectation that break-up only occurs in one direction whereas the solution remains regular in the other direction.
This is even more obvious in a close up in Fig. LABEL:derzoom (top), where we show in the region where the gradient catastrophe occurs. We recover the previously computed value in Fig. LABEL:derzoom (bottom), where is shown. It can be seen that , and we find that , as also predicted by the theoretical results in  for symmetry reasons.
The situation is similar for the dKP II case for the initial data (LABEL:uini1), which we treat analogously. We find that the fitting used for dKP I is optimal also for dKP II (), which yields the vanishing of at , , and . We observe that the solution develops a shock in the -direction () at time , see Fig. LABEL:u2ddkp2, where we show the solution to the dKP II equation with initial data (LABEL:uini1), at . Note that the tails with the algebraic fall off are now directed towards in contrast to the dKP I case. Thus the gradient catastrophe appears in this case first for negative (again we only treat the point of gradient catastrophe appearing first). By doing a fitting of , we find . The good agreement of this value with the maximum of the gradient of shows the self consistency of the approach.
In this case, the gradient catastrophe appears again only in one point , as we can see from Fig. LABEL:derivativesofu2kp2, where we show at , and .
The solutions to both dKP equations show thus a similar behavior, independent of the sign of in (LABEL:kpevol), except that of course, in the case of dKP I, the gradient catastrophe appears at , whereas in the case of dKP II, it appears at (as actually expected). As already mentioned, this is due to the tails with an algebraic fall off towards infinity which are for dKP I directed towards , whereas they are directed towards for dKP II.
3.5 KP solution in the small dispersion limit for initial data localized in one spatial dimension
Solutions to the KP equation will develop dispersive shocks, zones of rapid modulated oscillations, in the vicinity of a shock of solutions to the dKP equation for the same initial data. These were studied numerically for the first time in , see also . For KdV it was found that the difference between Hopf and KdV solution scales as for and as for . In this subsection, we establish numerically such scaling laws for the initial data localized in only one spatial direction.
For initial data of the form (LABEL:uini2), we show the solution to the KP I equation in the small dispersion limit with in Fig. LABEL:ut4s2d for several times. The computations are carried out with points for and .
As expected, the dispersive regularization of the break-up singularity leads to rapid modulated oscillations in the region where the corresponding dKP I solution develops a shock.
The number of oscillations increases as tends to , as already observed in . This can be seen in Fig. LABEL:ut4e0025 where the solution to the KP I equation in the small dispersion limit is shown for at .
In Fig. LABEL:ut4epss we present the solutions to the KP I equation in the small dispersion limit for several values of on the -axis for .
The contourplots of the solutions at are shown in Fig. LABEL:cont1 for several values of .
As usual, we ensure that the system is numerically well resolved by checking the decay of the Fourier coefficients, and the conservation of the numerically computed mass. The Fourier coefficients of the solution of the KP I equation in the small dispersion limit are plotted on the -axis for several values of in Fig. LABEL:coeftctmax at . They decrease to machine precision in all cases. The quantity is always smaller than , which indicates a numerical error well below plotting accuracy.
An important question is the scaling with of the norm of the difference between dKP and KP solutions for the same initial data. The norm of this difference is shown in Fig. LABEL:scaltc at (left) and at