## Abstract

The Cauchy problem for the Korteweg–de Vries (KdV) equation with small dispersion of order *ϵ*^{2}, *ϵ*≪1, is characterized by the appearance of a zone of rapid modulated oscillations. These oscillations are approximately described by the elliptic solution of KdV where the amplitude, wavenumber and frequency are not constant but evolve according to the Whitham equations. Whereas the difference between the KdV and the asymptotic solution decreases as *ϵ* in the interior of the Whitham oscillatory zone, it is known to be only of order *ϵ*^{1/3} near the leading edge of this zone. To obtain a more accurate description near the leading edge of the oscillatory zone, we present a multiscale expansion of the solution of KdV in terms of the Hastings–McLeod solution of the Painlevé-II equation. We show numerically that the resulting multiscale solution approximates the KdV solution, in the small dispersion limit, to the order *ϵ*^{2/3}.

## 1. Introduction

The mathematically rigorous study of the small dispersion limit of the Korteweg–de Vries (KdV) equation(1.1)with smooth initial data *u*_{0}(*x*) was initiated in the works of Lax & Levermore (1983), which stimulated intense research both numerically and analytically on the problem. The solution of the Cauchy problem of the KdV equation in the small dispersion limit is characterized by the appearance of a zone of rapid oscillations of frequency of order 1/*ϵ*, see for instance figure 1.

These oscillations are formed in the strong nonlinear regime and they have been analytically described in terms of elliptic functions, and in the general case in terms of theta functions in Venakides (1990) and Deift *et al*. (1997); the evolution in time of the oscillations was studied in Tian (1994). These results give a good asymptotic description of the oscillations only near the centre of the oscillatory zone (figure 2). In Grava & Klein (2007), we have studied numerically the small dispersion limit of the KdV equation for the concrete example of initial data of the form(1.2)We have compared the asymptotic description given in the works (Lax & Levermore 1983; Venakides 1990; Deift *et al*. 1997) with the numerical KdV solution. In Grava & Klein (2007), we have shown numerically that the difference between the KdV solution and the elliptic asymptotic solution at the centre of the oscillatory zone scales like *ϵ*, while this fails to be true at the boundary of the oscillatory zone. This fact was also observed for the Benjamin–Ono in Jorge *et al*. (1999). In particular at the left boundary, where the amplitude of the oscillations tends to zero, the difference between the KdV solution and the elliptic asymptotic solution scales like *ϵ*^{1/3}. In this manuscript we show that the Painlevé-II equation describes the envelope of the oscillations at the leading edge where the oscillations tend to zero. Painlevé equations appear in many branches of mathematics (for a review see Clarkson (2003)). For example, in the study of self-similar solutions of integrable equations, in the study of the Hele-Shaw flow near singularities (Fokas & Tanveer 1998) or in double-scaling limits in random matrix models (see for example Fokas *et al*. 1992; Brezin & Zee 1993; Bleher & Its 2003; Claeys *et al*. in press). In this work, following Kudashev & Suleimanov (1996), we perform a double-scaling limit of the KdV equation to derive the asymptotic description of the leading edge oscillations that are formed in the KdV small dispersion limit. We show that the envelope of the oscillations is determined by the Hastings & McLeod (1980) solution of the Painlevé-II equation. Then we compare numerically for the initial data *u*_{0}(*x*)=−sech^{2}*x* at the leading edge of the oscillatory front, the KdV solution with the derived multiscale solution and show that the difference between the two solutions scales like *ϵ*^{2/3}. We identify a neighbourhood of the leading edge of the Whitham zone where the multiscale solution gives a better asymptotic description than the asymptotic solution based on the elliptic and the Hopf solution. This allows the patching of different asymptotic descriptions to provide a more satisfactory treatment of the small dispersion limit of KdV as shown in figure 3.

Our analytical investigation of the multiscale expansion of the KdV solution requires the following assumptions on the initial data:

*u*_{0}(*x*) is negative with a single minimum; we chose the minimum to be at*x*=0 and normalized to −1, namely*u*_{0}(0)=−1;the function

*f*_{−}(*u*), which is the inverse of the monotone decreasing part of the initial data*u*_{0}(*x*), is such that ; and.

The last condition guarantees that the solution of the Cauchy problem for KdV exists for all times *t*>0.

This manuscript is organized as follows. In §2, we review the small dispersion limit of the KdV equation and study the small amplitude limit of the oscillations. In §3, we perform a multiscale expansion of the KdV solution when the oscillations tend to zero and show that the envelope of the oscillations is given by a solution of the Painlevé-II equation. In §4, we numerically compare the KdV solution with the multiscale solution obtained in §2. We show that the difference between the two solutions scales as *ϵ*^{2/3}, which is in accordance with our analytical result. We identify a zone near the leading edge where the multiscale solution provides a better description than the elliptic and the Hopf solution and patch the solutions. In §5, we summarize the results and add some concluding remarks on future directions of research.

## 2. Asymptotic solution of KdV in the small dispersion limit

The solution of the Cauchy problem for the KdV equation is characterized by the appearance of a zone of fast oscillations of wavelength of order *ϵ* (e.g. figure 1). These oscillations were called by Gurevich & Pitaevskii (1973) dispersive shock waves. Following the work of Lax & Levermore (1983), Venakides (1990) and Deift *et al*. (1997), the description of the small dispersion limit of the KdV equation is the following.

For 0≤

*t*<*t*_{c}, where*t*_{c}is a critical time, the solution*u*(*x*,*t*,*ϵ*) of the KdV Cauchy problem is approximated, for small*ϵ*, by*u*(*x*,*t*), which solves the Hopf equation(2.1)Here,*t*_{c}is the time when the first point of gradient catastrophe appears in the solution(2.2)of the Hopf equation. From the above, the time*t*_{c}of gradient catastrophe can be evaluated from the relationAfter the time of gradient catastrophe, the solution of the KdV equation is characterized by the appearance of an interval of rapid modulated oscillations. According to the Lax–Levermore theory, the interval [

*x*^{−}(*t*),*x*^{+}(*t*)] of the oscillatory zone is independent of*ϵ*. Here,*x*^{−}(*t*) and*x*^{+}(*t*) are determined from the initial data and satisfy the condition*x*^{−}(*t*_{c})=*x*^{+}(*t*_{c})=*x*_{c}, where*x*_{c}is the*x*-coordinate of the point of gradient catastrophe of the Hopf solution. Outside the interval [*x*^{−}(*t*),*x*^{+}(*t*)], the leading-order asymptotics of*u*(*x*,*t*,*ϵ*) as*ϵ*→0 is described by the solution of the Hopf equation (2.2). Inside the interval [*x*^{−}(*t*),*x*^{+}(*t*)], the solution*u*(*x*,*t*,*ϵ*) is approximately described, for small*ϵ*, by the elliptic solution of KdV (Gurevich & Pitaevskii 1973; Lax & Levermore 1983; Venakides 1990; Deift*et al*. 1997),(2.3)where(2.4)(2.5)with*K*(*s*) and*E*(*s*) the complete elliptic integrals of the first and second kind, ;*θ*is the Jacobi elliptic theta function defined by the Fourier series(2.6)

For constant values of the *β*_{i}, the formula (2.3) is an exact solution of KdV well known in the theory of finite-gap integration (Dubrovin & Novikov 1974; Its & Matveev 1975). However, in the description of the leading-order asymptotics of *u*(*x*, *t*, *ϵ*) as *ϵ*→0, the quantities *β*_{i} depend on *x* and *t* and evolve according to the Whitham equations (Whitham 1974)(2.7)where the speeds *v*_{i} are given by the formula(2.8)with *α* as in (2.5). The formula for *q* in the phase *Ω* in (2.4) that we give below was introduced in Grava & Klein (2007) and looks different but is equivalent to the one in Deift *et al*. (1997)(2.9)where *f*_{−}(*y*) is the inverse function of the decreasing part of the initial data *u*_{0}. The above formula for *q*(*β*_{1}, *β*_{2}, *β*_{3}) is valid as long as *β*_{1}>*β*_{2}>*β*_{3}>−1. When *β*_{3} reaches the minimum value 1 and passes over the negative hump, it is necessary to take into account also the increasing part of the initial data *f*_{+}(*u*) in formula (2.9). We denote by *T* this time. For *t*>*T*>*t*_{c} we introduce the variable *X*_{3} defined by *u*_{0}(*X*_{3})=*β*_{3} which is still monotonous. For values of *X*_{3} beyond the hump, namely *X*_{3}>0, we have to substitute (2.9) by the formula(2.10)The function *q*=*q*(*β*_{1}, *β*_{2}, *β*_{3}) is symmetric with respect to *β*_{1}, *β*_{2} and *β*_{3}, and satisfies a linear overdetermined system of Euler–Poisson–Darboux type. It has been introduced in the work of Tian (1994). The Whitham equations (2.7) can be integrated through the so-called hodograph transform that generalizes the method of characteristics and gives the solution in the implicit form (Tsarev 1985)(2.11)where the *v*_{i} are defined in (2.8) and the *w*_{i}=*w*_{i}(*β*_{1}, *β*_{2}, *β*_{3}) are obtained from an algebro-geometric procedure (Krichever 1988) by the formula (Tian 1994)(2.12)with *q* defined in (2.9) or (2.10). Formula (2.11) solves the initial-value problem for the Whitham equations (2.7) with boundary conditions

*leading edge*(2.13)*trailing edge*(2.14)

In Grava & Klein (2007), we have solved numerically the initial-value problem for the Whitham equations. In this way, we could perform a numerical comparison between the KdV small dispersion solution and the asymptotic formula (2.3) (figure 3). While in the interior of the oscillatory zone the error scales numerically like *ϵ*, at the left boundary of the oscillatory zone the error scales numerically like *ϵ*^{1/3}. To derive a more satisfactory asymptotic approximation of the KdV small dispersion limit in the vicinity of this point, we perform in §3 a double-scaling expansion of the KdV equation, following the double-scaling limits appearing in random matrix theory. Before doing this analysis, we study the elliptic solution (2.3) in the limit when the oscillations go to zero.

### (a) Small amplitude limit of the elliptic solution

We study the elliptic solution (2.3) near the leading edge, namely when oscillations go to zero. To avoid degeneracies, we rewrite the system (2.11) in the equivalent form(2.15)and perform the limit *δ*→0, whereTo simplify our calculation, we restrict ourselves to the case *t*_{c}<*t*<*T*. The following limit holds:(2.16)such that(2.17)Furthermore, the following identities hold:(2.18)(2.19)where(2.20)Substituting (2.17), (2.18) and (2.19) into (2.15), we arrive at the system(2.21)From the above, we deduce that, in the limit *δ*→0, the hodograph transform (2.15) reduces to the form (see Tian 1994; Grava & Tian 2002)(2.22)The above system enables one to determine *x*, *u* and *v* as a function of time. This time dependence will be denoted *x*=*x*^{−}(*t*), *u*=*u*(*t*) and *v*=*v*(*t*). We are interested in studying the behaviour of the elliptic solution (2.3) near the leading edge, namely when *x*=*x*^{−}(*t*) is small and *x*>*x*^{−}(*t*). For this purpose, we introduce two unknown functions of *x* and *t*,which tend to zero as *x*→*x*^{−}(*t*). We are going to derive the dependence of *Δ* as a function of *x*−*x*^{−}(*t*). Let us fix(2.23)Using the first equation of (2.21), we obtainExpanding the above expression near *β*_{1}(*x*, *t*)=*u*(*t*)+*Δ*(*x*, *t*), using the identity(2.24)and (2.22) we arrive at the expression(2.25)so that(2.26)Using the second equation in (2.21), we arrive at(2.27)where(2.28)Therefore,

*The elliptic solution* (*2.3*) *in the limit* (*2.23*) *takes the form*(2.29)*where the phase Ω*^{−} *takes the form*(2.30)*with*(2.31)*and u*(*t*)*, v*(*t*) *and x ^{−}*(

*t*)

*solve the system*(

*2.22*).

We first prove the relation (2.30). Using the expansionand (2.18), we obtain the following limit for the phase *Ω* in (2.4):Using the identity(2.22) and (2.19), we can rewrite the above in the formso that the phase *Ω* takes the form(2.32)We define(2.33)so that, expanding near *β*_{1}=*u*+*Δ*, by (2.22)(2.34)where *ϕ*_{2}(*x*, *t*) is defined in (2.31). To show that *η*_{0}(*u*, *v*) defined in (2.33) coincides with *ϕ*_{0} defined in (2.31), we differentiate (2.33) with respect to time,where we have used the identity (2.24) and . Integrating the r.h.s. of the above expression with respect to *t* from *t*_{c} to *t*, we obtain the formula (2.31).

Using (2.34) and (2.26), we rewrite the phase (2.32) in the form(2.35)where *ϕ*_{0} and *ϕ*_{2} are as in (2.31). Neglecting the higher order terms in *δ* and *Δ* of the above expansion one obtains (2.30).

Now we are ready to expand the theta-function expression in the limit of small amplitudes. Using (2.17) andone derives the small amplitude limit of the Jacobi *θ*-function (2.6)Substituting the above expansion in (2.3), one obtainswhich coincides with (2.29). ▪

## 3. Painlevé equations at the leading edge

In this section we present a multiscale description of the oscillatory behaviour of a solution to the KdV equation in the small dispersion limit close to the leading edge *x*^{−}(*t*) where *β*_{2}=*β*_{3}=*v* and *β*_{1}=*u*. We are interested in the double-scaling limit to the solution of the KdV equation (1.1) as *ϵ*→0 and *x*→*x*^{−}(*t*) in such a way that *x*−*x*^{−}(*t*)∝*ϵ*^{2/3}. We introduce the rescaled coordinate *y* near the leading edge,(3.1)which transforms the KdV equation (1.1), to the form(3.2)where . The substitution (3.1) has the effect that the linear term of (3.2) is just the Airy equation *u*_{t}+*u*_{yyy}=0 which has oscillatory solutions.

It is known (Kudashev & Suleimanov 1996; Deift *et al*. 1999) that the corrections to the Hopf solution near the leading edge are of the order *ϵ*^{1/3}. We thus make the ansatz(3.3)where *U*_{0}=*u*(*t*) is the solution at the leading edge. We assume that *U*_{k≥1} contains oscillatory terms with oscillations of the order 1/*ϵ*. In particular,(3.4)where(3.5)Similarly, we put(3.6)and(3.7)Terms proportional to sin(*ψ*/*ϵ*) can be absorbed by a redefinition of *ψ*. Since we impose no further restrictions on *ψ* here, such terms are therefore omitted in all orders. We consider only terms proportional to cos(*ψ*/*ϵ*) in order *ϵ*^{1/3} and the necessary terms in higher order to compensate the terms due to the nonlinearities in (3.2).

If we enter equation (1.1) with this ansatz, we immediately obtain from the term of order *ϵ*^{0} that *ψ*_{0,y}=*ψ*_{1,y}=0. From the term of order *ϵ*^{1/3} we get(3.8)In order *ϵ*^{2/3} we obtain the following equations:(3.9)(3.10)(3.11)In order *ϵ* we get(3.12)(3.13)(3.14)(3.15)(3.16)(3.17)A solution to (3.10), (3.11) and (3.14) is(3.18)where *C*(*t*) is a free function of *t*. Note that(3.19)where *u*(*t*) and *v*(*t*) are defined in (2.22). The above equations imply thatwhere *h*(*t*) is a free function of *t*. Comparing the above relation with the formula of *ϕ*_{2} in (2.29) of the phase in the small amplitude expansion, we can conclude that *C*(*t*)=0 and *h*(*t*)=0 so that(3.20)From (3.8), (3.19) and (3.20), we derive thatnamely(3.21)From (3.12) we find(3.22)where *k*(*t*) is a free function of *t*. It will be fixed by matching it with the elliptic solution in the Whitham zone. Substituting (3.22) and (3.20) into (3.15), we arrive at the equation(3.23)Making the substitution(3.24)with *y*_{0}=12*k*(*t*)(*u*(*t*)−*v*(*t*))/*v*_{t}, we arrive at the equation(3.25)which is a special case of the Painlevé-II equation *A*_{zz}=*zA*+2*A*^{3}−*γ*, with *γ*=0.

Since we are interested only in terms up to order *ϵ*^{1/3} in *u*(*x*, *t*, *ϵ*), the terms *b*_{1}, *b*_{2}, *c*_{0}, *c*_{2} and *c*_{3} are not important for us. However, we had to go to order *ϵ* to determine *ψ*_{3}, which will contribute to the *ϵ*^{1/3} terms in *u*.

To sum up, we get for *u*(*x*, *t*, *ϵ*)(3.26)whereand *h* is an integration constant. There are free functions of *t* in the integration of the multiscale equations, namely the functions *k*(*t*) and *ψ*_{3}(*t*) and the constant *h*. Moreover, the solution of the Painlevé-II equations needs to be fixed. We fix the constants by comparing the multiscale solution (3.28) with the elliptic solution (2.29) at the border of the Whitham zone. Indeed, comparing (3.26), (2.29) and (2.35) we obtain(3.27)and *k*(*t*)=0, *ψ*_{3}(*t*)=0 and *h*=0. Therefore, the multiscale solution takes the form(3.28)where *y*=*ϵ*^{−2/3}(*x*−*x*^{−}(*t*)), *a*(*y*, *t*) satisfies (3.23) andFor the numerical comparison in the following section, we consider terms up to order *ϵ*^{1/3} in *u* in (3.28):(3.29)where *ψ*(*y*, *t*) is as given above. For fixing the particular solution of the Painlevé-II equation (3.25), the following considerations are needed. For large *x*<*x*^{−}(*t*), the solution of KdV is essentially approximated by the solution of the Hopf equation, and the term of order *ϵ*^{1/3} has to be negligible in (3.29), namely for large negative *y*=(*x*−*x*^{−}(*t*))*ϵ*^{−2/3}. For *x*<*x*^{−}(*t*),because *v*_{t}=6/((*u*−*v*)∂_{vv}*Φ*(*v*;*u*))<0 since ∂_{vv}*Φ*(*v*;*u*)<0 and *u*>*v*. Therefore, for *x*≪*x*^{−}(*t*), we conclude that(3.30)For *x*>*x*^{−}(*t*), from the small amplitude limit of the elliptic solution of the KdV equation, we obtain by combining (2.27), (2.28) and (3.27)which, in the limit *ϵ*→0 or *y*→+∞, givesor equivalently, by (3.24),(3.31)The existence and uniqueness of the solution of (3.25) satisfying (3.30) and (3.31) was first established by Hastings & McLeod (1980; see also later works of Kapaev 1992; Clarkson 2003). It is also worth noticing that the asymptotics of *A*(*z*) at *z*→+∞can be specified as(3.32)where Ai(*z*) is the Airy function. Moreover, the asymptotic condition (3.32) characterizes the solution *A*(*z*) uniquely, so that (3.31) and (3.32) constitute an example of the so-called connection formula for the Painlevé equations (e.g. Fokas & Its 1993; Clarkson & McLeod 1988).

## 4. Comparison of the multiscale expansion and the asymptotic solution to the small dispersion KdV

The numerical evaluation of the asymptotic solution based on the Hopf and the elliptic solution is described in Grava & Klein (2007). To evaluate the multiscale solution (3.28), one needs, in addition to the quantities computed there, the Hastings–McLeod solution to the Painlevé-II equation. This solution was calculated numerically by Tracy & Widom (1994) with standard solvers for the ordinary differential equation and by Prähofer & Spohn (2004; www-m5.ma.tum.de/KPZ/) with in principle arbitrary precision with a Taylor series approach. The general family of solutions of Painlevé-II such that with *a* positive constant was studied numerically in Rosales (1978) and analytically in Clarkson & McLeod (1988). Solutions to Painlevé-II in the complex plane were studied analytically and numerically by Fokas & Tanveer (1998). We use here an approach based on spectral methods which is described briefly in appendix A. This approach is both efficient and of high precision and can directly be combined with the numerics of Grava & Klein (2007).

### (a) Times *t*≫*t*_{c}

Close to break-up the multiscale expansion is expected to be inefficient since it is best near the leading edge, and since at break-up both the leading and the trailing edge coincide. We will discuss this solution close to break-up below, but first we will study it for time *t*=0.4≫*t*_{c}=0.216…. In figure 4, one can see that the multiscale solution gives an excellent approximation of the KdV solution for *x*<*x*^{−}(0.4)=−3.2297 and in the Whitham zone close to *x*^{−}. For larger values of *x*, the solutions are out of phase and the values of the multiscale solution are shifted towards positive values.

The difference of the two solutions is shown in figure 5. From this figure it is even more obvious that the multiscale solution is a valid approximation in the Whitham zone near the leading edge, but the difference increases rapidly for |*x*|≫*x*^{−}.

### (b) *ϵ* dependence

In Grava & Klein (2007), it was shown that the asymptotic description becomes more accurate with decreasing *ϵ*. The same is true for the multiscale solution as can be seen in figure 6. The zone, in which the multiscale solution gives a better approximation than the asymptotic elliptic solution, shrinks with *ϵ*. For *x*≫*x*^{−}(*t*), the multiscale solution is always only a poor approximation to the KdV solution. The maximal difference *Δ*_{max} of the KdV solution and the multiscale solution near this edge decreases roughly as *ϵ*^{2/3}. More precisely, the error can be fitted with a straight line by a standard linear regression analysis, −log *Δ*_{max}=−*a* log *Δ*_{max}=−*a* log *ϵ*+*b* with *a*=0.63, *b*=0.41. The correlation coefficient is *r*=0.999, the standard error is *σ*_{a}=0.02.

### (c) Comparison and matching with the asymptotic solution

The aim of this paper is to improve the asymptotic description of the small dispersion limit of KdV near the leading edge. In figure 7, it can be seen that the multiscale solution will indeed be a much better approximation near this edge. Near the leading edge, the multiscale solution provides a superior description of the KdV solution, whereas the elliptic asymptotic solution is much better for *x*≫*x*^{−}(*t*) in the Whitham zone. In fact it is possible to identify a zone where the multiscale solution is more satisfactory than the asymptotic solution. Owing to the strong oscillations of the solutions, there is a certain ambiguity in the definition of this zone. We define the limits of the zone as the last intersection (or where the solutions come closest) on which the other solution has an error with larger oscillations. In this zone it is possible to replace the asymptotic solution by the multiscale solution. The result of this patchwork approach is shown in figure 8. It can be seen that the resulting amended asymptotic description has an accuracy near the leading edge of the same order as in the interior of the Whitham zone. The maximal difference between the KdV and the asymptotic solution still occurs near the leading edge.

As already mentioned, the zone where the multiscale solution provides a better approximation to the KdV solution than the asymptotic solution shrinks with *ϵ* as can be inferred from figure 9. The width of this zone decreases roughly as *ϵ*^{2/3}, which shows the self-consistency of the used rescaling of the spatial coordinate near the leading edge. More precisely, we find a scaling *ϵ*^{a} with *a*=0.66, correlation coefficient *r*=0.9996 and standard error *σ*_{a}*=*0.015. It can be seen that the zone is not symmetric around the leading edge; it extends much further into the Hopf region than in the Whitham zone. This is due to the fact that the multiscale solution is quickly out of phase with the rapid oscillations in the Whitham zone, and that the Hopf solution does not have oscillations.

### (d) Break-up time

In Grava & Klein (2007), it was shown that the elliptic asymptotic solution is worst near the break-up of the Hopf solution. The multiscale expansion obtained in the previous section is not defined for times before *t*_{c}, and it will be worst there since it can be understood as an expansion around the leading edge of the Whitham zone. At break-up, however, leading and trailing edges coincide. Thus, the approximation is rather crude there, but it increases in quality with time as can be seen in figure 10. It is, however, interesting to study at which times the multiscale solution starts to give a better asymptotic description than other approaches. Dubrovin conjectured (Dubrovin 2006) that the asymptotic behaviour of the KdV solution close to the break-up of the corresponding Hopf solution is given by a particular solution to the second equation in the Painlevé-I hierarchy. In Grava & Klein (2006), we provided strong numerical evidence for the validity of this conjecture. The natural question is whether Painlevé-II description near the critical point provides a satisfactory asymptotic solution for KdV until times where the multiscale solution studied in the present paper provides a valid description near the leading edge. A comparison of figure 10 with a similar figure in Grava & Klein (2006) shows that this is qualitatively the case. In figure 10 the multiscale solution is shown for *x*<*x*^{+}(*t*). Near break-up the approximation is acceptable only close to the break-up point. For larger times, more and more oscillations are satisfactorily reproduced by the multiscale solution. But the multiscale solution will be a better approximation of the oscillations than the Hopf solution only for times *t*≫*t*_{c}.

## 5. Outlook

In the present work we have considered a multiscale solution to the KdV equation in the small dispersion limit close to the leading edge of the oscillatory zone. We studied the solution up to order *ϵ*^{1/3} and show that the envelope of the oscillations is described by the Hastings–McLeod solution of the Painlevé-II. The validity of the approach in the considered limit was shown numerically. The double-scaling expansion of the KdV solution in the small dispersion limit will be investigated with the Riemann–Hilbert approach and the steepest descent method for the oscillatory Riemann–Hilbert problem as done in Deift *et al*. (1997). The Riemann–Hilbert approach seems so far the only analytical tool to study the double-scaling expansion to the *Cauchy problem* of the KdV equation. This project will be the subject of our future research.

As can be seen from figure 3, the asymptotic solution of the KdV equation does not give a satisfactory description of the KdV small dispersion limit also at the trailing edge of the oscillatory zone. This problem will be investigated in a subsequent publication, too. A similar problem was tackled in the context of matrix models in Claeys (submitted).

## Acknowledgments

We thank B. Dubrovin and J. Frauendiener for their helpful discussions and hints. We acknowledge the support by the MISGAM programme of the European Science Foundation. T.G. acknowledges support by the RTN ENIGMA and Italian COFIN 2004 ‘Geometric methods in the theory of nonlinear waves and their applications’. The authors wish to thank the referees for the improvements suggested to the manuscript.

## Numerical solution of the Painlevé-II equation

We are interested in the numerical computation of the Hastings–McLeod solution to the Painlevé-II equation(A1)which is subject to the asymptotic conditions (Hastings & McLeod 1980)(A2)(A3)where Ai(*z*) is the Airy function. Numerically, we will consider equation (A 1) on a finite interval [*z*_{l},*z*_{r}] (typically [−10,10]). The asymptotic solution near ±∞, which will be discussed in more detail below, is truncated in a way that the truncation error at *z*_{l}, *z*_{r} is below 10^{−10}. At these points, we impose the values following from the asymptotic solutions as boundary conditions, namely

(A4)

To solve equation (A 1) for *z*∈[*z*_{l},*z*_{r}], we use spectral methods since they allow for an efficient numerical approximation of high accuracy. We map the interval [*z*_{l},*z*_{r}] with a linear transformation *z*→*x* to the interval *I*=[−1,1] and expand *A* there in Chebyshev polynomials.

Let us briefly summarize the Chebyshev approach (for details see, for example, Canuto *et al*. 1988; Fornberg 1996; Frauendiener & Klein 2004). The Chebyshev polynomials *T*_{n}(*x*) are defined on the interval *I* by the relationA function *f* on *I* is approximated via Chebyshev polynomials, where the spectral coefficients *a*_{n} are obtained by the conditions , *l*=0, …, *N*. This approach is called a collocation method. If the collocation points are chosen to be *x*_{l}=cos(*πl*/*N*), the spectral coefficients follow from *f* via a discrete cosine transform (DCT) for which fast algorithms exist. We use here a DCT within Matlab. A recursive relation for the derivative of Chebyshev polynomials implies that the action of the differential operator ∂_{x} on *f*(*x*) leads to an action of a matrix *D* on the vector of the spectral coefficients *a*_{n}. Thus, we express *A*(*x*) in terms of Chebyshev polynomials, (we typically work with *N*=128), and the coefficients of ∂_{x}*A* in terms of Chebyshev polynomials are determined then via .

To solve equation (A 1) on the interval [*z*_{l},*z*_{r}], we use an iterative approach,(A5)We start with . In each step of the iteration, we solve equation (A 5) for *A*_{n+1} with the boundary conditions (A 4). The boundary conditions are imposed with a *τ*-method: the last two rows of the matrix *D*^{2} for the second derivative are replaced with the boundary conditions at *x*=±1. Since *T*_{n}(±1)=(±1)^{n}, the resulting matrix *L*, which will be inverted in each step of the iteration, has only 1 and −1 as entries in the last two rows and is thus better conditioned than the matrix *D*^{2}. It turns out that the iteration is unstable if no relaxation is used. We thus define with *μ*=0.009. With this choice of the parameters, the iteration converges. It is stopped when the difference between *A*_{n+1} and *A*_{n} is of the order of machine precision (Matlab works internally with a precision of the order of 10^{−16}; owing to rounding errors, machine precision is typically limited to the order of 10^{−14}). The solution is shown in figure 11.

To test the accuracy of the solution, we plot in figure 12 the quantity *P*_{II}*A* as computed with spectral methods on the collocation points. It can be seen that the error is the biggest on the boundary which is even more obvious from figure 13. The found solution is also compared with a numerical solution with a standard ode solver as *bvp4c* in Matlab. The solutions agree within the limits of numerical precision.

For general values of *z*, the solution is obtained as follows: for values of *z*∈[*z*_{l},*z*_{r}] they follow from the spectral data via . Note that the accuracy of the solution is the best on the collocation points, but we can expect it to be of the order of at least 10^{−6} even at points *z* in between. For values of *z*<*z _{l}*, we use the approximation ; for values of

*z*>

*z*

_{r}, we use the approximation . This provides a global approximation to the solution with an accuracy of the order of 10

^{−6}and better, which is sufficient for our purposes. Higher precision can be reached within the used approach without problems: one can either increase the values of −

*z*

_{l}and

*z*

_{r}and use a higher number of polynomials, or use higher order terms in the asymptotic solution of

*A*for

*z*→±∞.

## Footnotes

- Received June 19, 2007.
- Accepted November 29, 2007.

- © 2008 The Royal Society