---
abstract: 'We study numerically solutions to the Korteweg-de Vries and Camassa-Holm equation close to the breakup of the corresponding solution to the dispersionless equation. The solutions are compared with the properly rescaled numerical solution to a fourth order ordinary differential equation, the second member of the Painlevé I hierarchy. It is shown that this solution gives a valid asymptotic description of the solutions close to breakup. We present a detailed analysis of the situation and compare the Korteweg-de Vries solution quantitatively with asymptotic solutions obtained via the solution of the Hopf and the Whitham equations. We give a qualitative analysis for the Camassa-Holm equation'
address:
- 'SISSA, via Beirut 2-4, 34014 Trieste, Italy'
- 'Max Planck Institute for Mathematics in the Sciences, Inselstr. 22, 04103 Leipzig, Germany'
author:
- Tamara Grava
- Christian Klein
title: ' Numerical study of a multiscale expansion of KdV and Camassa-Holm equation'
---

[^1]

Introduction
============

It is well known that the solution of the Cauchy problem for the Hopf equation $$\label{hopf}
u_t+6uu_x=0,\quad u(x,t=0)=u_0(x), \;\;x\in\mathbb{R},\;t\in\mathbb{R}^+$$ reaches a point of gradient catastrophe in a finite time. The solution of the viscosity or conservative regularization of the above hyperbolic equation display a considerably different behavior. Equation (\[hopf\]) admits an Hamiltonian structure $$u_t +\{ u(x), H_0\}\equiv u_t +\partial_x\frac{\delta H_0}{\delta u(x)} =0,\quad$$ with Hamiltonian and Poisson bracket $$H_0 =\int u^3\, dx, \quad
\{ u(x) , u(y)\}=\delta'(x-y),$$ respectively. All the Hamiltonian perturbations up to the order $\epsilon^4$ of the hyperbolic equation (\[hopf\]) have been classified in [@dubcr]. They are parametrized by two arbitrary functions $c(u)$, $p(u)$ $$\begin{split}
\label{riem2}
&
u_t +6u\, u_x + \frac{\epsilon^2}{24} \left[ 2 c\, u_{xxx} + 4 c' u_x u_{xx}
+ c'' u_x^3\right]+\epsilon^4 \left[ 2 p\, u_{xxxxx} \right.\\
&
\\
&\left.
+2 p'( 5 u_{xx} u_{xxx} + 3 u_x u_{xxxx}) + p''( 7 u_x u_{xx}^2 + 6 u_x^2 u_{xxx} ) +2 p''' u_x^3 u_{xx}\right]=0,
\end{split}$$ where the prime denotes the derivative with respect to $u$. The corresponding Hamiltonian takes the form $$H=\int \left[ u^3 - \epsilon^2 \frac{c(u)}{24} u_x^2
+\epsilon^4  p(u) u_{xx}^2 \right]\, dx$$ For $c(u)=12$, $p(u)=0$ one obtains the Korteweg - de Vries (KdV) equation $u_t+6uu_x+\epsilon^{2}u_{xxx}=0$, and for $c(u)=48u $ and $p(u)=2u$ the Camassa-Holm equation up to order $\epsilon^{4}$; for generic choices of the functions $c(u)$, $p(u)$ equation (\[riem2\]) is apparently not an integrable PDE. However it admits an infinite family of commuting Hamiltonians up to order $O({\epsilon}^6).$

The case of small viscosity perturbations of one-component hyperbolic equations has been well studied and understood (see [@bressan] and references therein), while the behavior of solutions to the conservative perturbation (\[riem2\]) to the best of our knowledge has not been investigated after the point of gradient catastrophe of the unperturbed equation except for the KdV case, [@LL; @V2; @DVZ].

In a previous paper [@GK] (henceforth referred to as I) we have presented a quantitative numerical comparison of the solution of the Cauchy problem for KdV $$\label{KdV}
u_t+6uu_x+\epsilon^{2}u_{xxx}=0,\quad u(x,0)=u_0(x),$$ in the small dispersion limit $\epsilon{\rightarrow}0$, and the asymptotic formula obtained in the works of Lax and Levermore [@LL], Venakides [@V2] and Deift, Venakides and Zhou [@DVZ] which describes the solution of the above Cauchy problem at the leading order as ${\epsilon}\rightarrow 0$. The asymptotic description of [@LL],[@DVZ] gives in general a good approximation of the KdV solution, but is less satisfactory near the point of gradient catastrophe of the hyperbolic equation. This problem has been addressed by Dubrovin in [@dubcr], where, following the universality results obtained in the context of matrix models by Deift et all [@DKMVZ], he formulated the universality conjecture about the behavior of a generic solution to the Hamiltonian perturbation (\[riem2\]) of the hyperbolic equation (\[hopf\]) near the point $(x_c,t_c,u_c)$ of gradient catastrophe for the solution of (\[hopf\]). He argued that, up to shifts, Galilean transformations and rescalings, this behavior essentially depends neither on the choice of solution nor on the choice of the equation. Moreover, the solution near the point $(x_c, t_c, u_c)$ is given by $$\label{univer}
u(x,t,{\epsilon})\simeq u_c +a\,\epsilon^{2/7} U \left( b\, \epsilon^{-6/7} (x- x_c-6u_c (t-t_c)); c\, 
\epsilon^{-4/7} (t-t_c)\right) +O\left( \epsilon^{4/7}\right)$$ where $a$, $b$, $c$ are some constants that depend on the choice of the equation and the solution and $U=U(X; T)$ is the unique real smooth solution to the fourth order ODE $$\label{PI2}
X=6T\, U -\left[ U^3 + (\frac12 U_{X}^2 +  U\, U_{XX} ) +\frac1{10} U_{XXXX}\right],$$ which is the second member of the Painlevé I hierarchy. We will call this equation PI2. The relevant solution is characterized by the asymptotic behavior $$\label{PI2asym}
        U(X,T)=\mp(X)^{\frac{1}{3}}\mp \dfrac{2T}{X^{\frac13}}+O(X^{-1}),\quad X{\rightarrow}\pm \infty,$$ for each fixed $T\in\mathbb{R}$. The existence of a smooth solution of (\[PI2\]) for all $X,T\in\mathbb{R}$ satisfying (\[PI2asym\]) has been recently proved by Claeys and Vanlessen [@CV]. Furthermore they study in [@CV1] the double scaling limit for the matrix model with the multicritical index and showed that the limiting eigenvalues correlation kernel is obtained from the particular solution of (\[PI2\]) satisfying (\[PI2asym\]). This result was conjectured in the work of Brézin, Marinari and Parisi [@BMP].

In this paper we address numerically the validity of (\[univer\]) for the KdV equation, and we identify the region where this solution provides a better description than the Lax-Levermore, and Deift-Venakides-Zhou theory. As an outlook for the validity of (\[univer\]) for other equations in the family (\[riem2\]), we present a numerical analysis of the Camassa-Holm equation near the breakup point. While the validity of (\[univer\]) can be theoretically proved using a Riemann-Hilbert approach to the small dispersion limit of the KdV equation [@DVZ] and recent results in [@DKMVZ],[@CV],[@CV1], for the Camassa-Holm equation and also for the general Hamiltonian perturbation to the hyperbolic equation (\[hopf\]), the problem is completely open. Furthermore for the general equation (\[riem2\]), the existence of a smooth solution for a short time has not been established yet. An equivalent analysis should also be performed for Hamiltonian perturbation of elliptic systems, in particular for the semiclassical limit of the focusing nonlinear Schrödinger equation [@KMM],[@TVZ].

The paper is organized as follows. In section 2 we give a brief summary of the Lax-Levermore, and Deift-Venakides-Zhou theory and the multiscale expansion (\[univer\]). In section 3 we present the numerical comparison between the asymptotic description based on the Hopf and Whitham solutions and the multiscale solutions with the KdV solution. In section 4 we study the same situation for the Camassa-Holm equation. In the appendix we briefly outline the used numerical approaches.

Asymptotic and multiscale solutions
===================================

Following the work of [@LL], [@V2] and [@DVZ], the rigorous theoretical description of the small dispersion limit of the KdV equation is the following: Let $\bar{u}(x,t)$ be the zero dispersion limit of $u(x,t,\epsilon)$, namely $$\label{baru}
\bar{u}(x,t)=\lim_{{\epsilon}{\rightarrow}0}u(x,t,{\epsilon}).$$ 1) for $0\leq t< t_c$, where $t_c$ is a critical time, the solution $u(x,t,\epsilon)$ of the KdV Cauchy problem is approximated, for small ${\epsilon}$, by the limit $\bar{u}(x,t)$ which solves the Hopf equation $$\label{Hopf}
\bar{u}_t+6\bar{u}\bar{u}_x=0.$$ Here $t_c$ is the time when the first point of gradient catastrophe appears in the solution $$\label{Hopfsol}
\bar{u}(x,t)=u_0(\xi),\quad x=6tu_0(\xi)+\xi,$$ of the Hopf equation. From the above, the time $t_c$ of gradient catastrophe can be evaluated from the relation $$t_{c}=\dfrac{1}{\min_{\xi\in\mathbb{R}}[-6u_0'(\xi)]}.$$ 2) After 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 $\epsilon$. 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,{\epsilon})$ as ${\epsilon}{\rightarrow}0$ is described by the solution of the Hopf equation (\[Hopfsol\]). Inside the interval $[x^-(t), x^+(t)]$ the solution $u(x,t,{\epsilon})$ is approximately described, for small ${\epsilon}$, by the elliptic solution of KdV [@GP], [@LL], [@V2], [@DVZ], $$\label{elliptic}
u(x,t,{\epsilon})\simeq \bar{u}+ 2{\epsilon}^2\frac{\partial^2}{\partial
x^2}\log\theta\left(\dfrac{\sqrt{\beta_1-\beta_3}}{2{\epsilon}K(s)}[x-2 t(\beta_1+\beta_2+\beta_3) -q];\mathcal{T}\right)$$ where now $\bar{u}=\bar{u}(x,t)$ takes the form $$\label{ubar}
\bar{u}=\beta_1+\beta_2+\beta_3+2\alpha,$$ $$\label{alpha}
\alpha=-\beta_{1}+(\beta_{1}-\beta_{3})\frac{E(s)}{K(s)},\;\;\mathcal{T}=i\dfrac{K'(s)}{K(s)},
\;\; s^{2}=\frac{\beta_{2}-\beta_{3}}{\beta_{1}-\beta_{3}}$$ with $K(s)$ and $E(s)$ the complete elliptic integrals of the first and second kind, $K'(s)=K(\sqrt{1-s^{2}})$; $\theta$ is the Jacobi elliptic theta function defined by the Fourier series $$\theta(z;\mathcal{T})=\sum_{n\in\mathbb{Z}}e^{\pi i n^2\mathcal{T}+2\pi i nz}.$$ For constant values of the $\beta_i$ the formula (\[elliptic\]) is an exact solution of KdV well known in the theory of finite gap integration [@IM], [@DN0]. However in the description of the leading order asymptotics of $u(x,t,{\epsilon})$ as ${\epsilon}{\rightarrow}0$, the quantities $\beta_i$ depend on $x$ and $t$ and evolve according to the Whitham equations [@W] $$\dfrac{\partial}{\partial t}\beta_i+v_i\dfrac{\partial}{\partial x}\beta_i=0,\quad i=1,2,3,$$ where the speeds $v_i$ are given by the formula $$v_{i}=4\frac{\prod_{k\neq
     i}^{}(\beta_{i}-\beta_{k})}{\beta_{i}+\alpha}+2(\beta_1+\beta_{2}+\beta_{3}),
    \label{eq:la0}$$ with $\alpha$ as in (\[alpha\]). Lax and Levermore first derived, in the oscillatory zone, the expression (\[ubar\]) for $\bar{u}=\bar{u}(x,t)$ which clearly does not satisfy the Hopf equation. The theta function formula (\[elliptic\]) for the leading order asymptotics of $u(x,t,{\epsilon})$ as ${\epsilon}{\rightarrow}0$, was obtained in the work of Venakides and the phase $q=q(\beta_1,\beta_2,\beta_3)$ was derived in the work of Deift, Venakides and Zhou [@DVZ], using the steepest descent method for oscillatory Riemann-Hilbert problems [@DZh] $$\label{q0}
    q(\beta_{1},\beta_{2},\beta_{3}) = \frac{1}{2\sqrt{2}\pi}
    \int_{-1}^{1}\int_{-1}^{1}d\mu d\nu \frac {f_-( \frac{1+\mu}{2}(\frac{1+\nu}{2}\beta_{1}
    +\frac{1-\nu}{2}\beta_{2})+\frac{1-\mu}{2}\beta_{3})}{\sqrt{1-\mu}
    \sqrt{1-\nu^{2}}},$$ where $f_-(y)$ is the inverse function of the decreasing part of the initial data. The above formula holds till some time $T>t_c$ (see [@DVZ] or I for times $t>T$).

3\) Fei-Ran Tian proved that the description in 1) and 2) is generic for some time after the time $t_c$ of gradient catastrophe [@FRT1].

In I we discussed the case $u_{0}(x) = -\mbox{sech}^{2}x$ in detail as an example. The main results were that the asymptotic description is of the order $\mathcal{O}(\epsilon)$ close to the center of the Whitham zone, but that the approach gives considerably less satisfactory results near the edges of the Whitham zone and close to the breakup of the corresponding solution to the Hopf equation. In the present paper we address the behavior near the point of gradient catastrophe of the Hopf solution in more detail. In Fig. \[fig2\] we show the KdV solution and the corresponding asymptotic solution as given above for several values of the time near the critical time $t_{c}$. It can be seen that there are oscillations before $t_{c}$, and that the solution in the Whitham zone provides only a crude approximation of the KdV solution for small $t-t_{c}$.

The situation does not change in principle if we consider smaller values of $\epsilon$ as can be seen from Fig. \[figbreak4\]. The solution shows the same qualitative behavior as in Fig. \[fig2\], just on smaller scales in $t$ and $x$.

Multiscale expansion
--------------------

We give a brief summary of the results in [@dubcr] relevant for the KdV case we are interested in here. Near the point of gradient catastrophe $(x_c,t_c, u_c)$, the Hopf solution is generically given in lowest order by the cubic $$\label{hopfsolc}
x-x_c \simeq 6(t-t_c) u -k(u-u_c)^3,\quad k=-f_-'''(u_c)/6,$$ because $6t_c+f_-'(u_c)=0$ and $f_c''(u_c)=0$. Here $f_-(u)$ is the inverse of the decreasing part of the initial data $u_0(x)$. Now let us consider $h_k=\dfrac{\delta H_k}{\delta u}$ where $H_k$ are the KdV Hamiltonians such that $h_k=u^{k+2}/(k+2)!+O(\epsilon^2)$. We have $$h_{-1}=u,\;\;h_0=\dfrac{u^2}{2}+\dfrac{\epsilon^2}{6} u_{xx},\;\; h_1=\dfrac{1}{6}(u^3+\frac{\epsilon^2}{2}(u_x^2+2uu_{xx})+\dfrac{\epsilon^4}{10}u_{xxxx}),$$ and the KdV equation is obtained from $u_t+6\partial_xh_0=0$. Then $$x=6tu+a_0h_0+a_1h_1+\dots a_kh_k,$$ is a symmetry of the KdV equation [@DZ]. Setting $a_0=0,$ $a_1=-f'''(u_c)/6=k$ and $a_{k>2}=0$, and making the shift $t\rightarrow t-t_c$, $u\rightarrow u-u_c$ and the Galilean transformation $x\rightarrow x-x_c-6(t-t_c)u_c$ we arrive at the fourth order equation of Painlevé type $$\label{painleve}
x-x_c-6(t-t_c)u_c = 6(t-t_c) (u-u_c) -k\left[(u-u_c)^3 +\epsilon^2
(\dfrac{ u_x^2}{2} +(u-u_c) u_{xx})+\dfrac{ \epsilon^4 }{10} u_{xxxx}\right]$$ which is an exact solution of the KdV equation and can be considered as a perturbation of the Hopf solution (\[hopfsolc\]) near the point of gradient catastrophe $(x_c,t_c,u_c)$. The solution $u(x,t,{\epsilon})$ of (\[painleve\]) is related to the solution $U(X,T)$ of (\[PI2\]) by the rescalings $$\label{KdVrescaled}
u(x,t,{\epsilon})=u_c+\left(\dfrac{{\epsilon}}{  k}\right)^{2/7} U(X,T)$$ where $$\label{scalings}
X =\dfrac{ x - x_c - 6 u_c (t-t_c)}{\epsilon^{\frac{6}{7}}k^{\frac17}},\quad 
T=\dfrac{t-t_c}{{\epsilon}^{\frac{4}{7}}k^{\frac37}}.$$ According to the conjecture in [@dubcr], the solution (\[KdVrescaled\]) is an approximation modulo terms $O({\epsilon}^{\frac{4}{7}})$ to the solution of the Cauchy problem (\[KdV\]) for $(x,t,u)$ near the point of gradient catastrophe $(x_c,t_c,u_c)$ of the hyperbolic equation (\[Hopf\]).

Numerical comparison
====================

In this section we will present a comparison of numerical solutions to the KdV equation and asymptotic solutions arising from solutions to the Hopf and the Whitham equations as well as the Painlevé I2 equation as given above. Since we control the accuracy of the used numerical solutions, see I, [@numart1d] and the appendix, we ensure that the presented differences are entirely due to the analytical description and not due to numerical artifacts. We study the $\epsilon$-dependence of these differences by linear regression analysis. This will be done for nine values of $\epsilon$ between $10^{-1}$ and $10^{-3}$. Obviously the numerical results are only valid for this range of parameters, but it is interesting to note the high statistical correlation of the scalings we observe. We consider the initial data $$u_0(x)=-1/\cosh^2x.$$ For this initial data $$x_c=-\dfrac{\sqrt{3}}{2}+\log((\sqrt{3}-1)/\sqrt{2}),
\;\;t_c=\dfrac{\sqrt{3}}{8},\;\;
   \;\;u_c=-2/3.
    \label{tcrit}$$

Hopf solution
-------------

We will first check whether the rescalings of the coordinates given in (\[KdVrescaled\]) are consistent with the numerical results. It is known that the Hopf solution provides for times $t\ll t_{c}$ an asymptotic description of the KdV solution up to an error of the order $\epsilon^{2}$. This means that the $L_{\infty}$-norm of the difference between the two solutions decreases as $\epsilon^{2}$ for $\epsilon{\rightarrow}0$. For $t=0.1$ we actually observe this dependence. More precisely this difference $\Delta_{\infty}$ can be fitted with a straight line by a standard linear regression analysis, $-\log_{10}\Delta_{\infty}=-a\log_{10}\epsilon+b$ with $a=1.9979$, with a correlation coefficient of $r=0.99999$ and standard error $\sigma_{a}=4.1*10^{-3}$.

Near the critical time $t_{c}$ this picture is known to change considerably. Dubrovin’s conjecture [@dubcr] presented above suggests that the difference between Hopf and KdV solution near the critical point should scale roughly as $\epsilon^{2/7}$. In the following we will always compare solutions in the intervals $$\label{interval}
[x_{c}+6u_{c}(t-t_{c})-\alpha 
\epsilon^{6/7},x_{c}+6u_{c}(t-t_{c})+\alpha \epsilon^{6/7}]$$ where $\alpha$ is an $\epsilon$-independent constant (typically we take $\alpha=3$).

Numerically we find at the critical time that the $L_{\infty}$-norm of the difference between Hopf and KdV solution scales like ${\epsilon}^a$ where $a=0.2869$ ($2/7=0.2857\ldots$) with correlation coefficient $r=0.9995$ and standard error $\sigma_{a}=6.9*10^{-3}$. Thus we confirm the expected scaling behavior within numerical accuracy. We also test this difference for times close to $t_{c}$. The relations (\[KdVrescaled\]) suggest, however, a rescaling of the time, i.e., to compare solutions for different values of $\epsilon$ at the same value of $T$. We compute the respective solutions for KdV times $t_{\pm}(\epsilon)=t_{c}\pm 0.1\epsilon^{4/7}$. Before breakup at $t_{-}$ we obtain $a=0.31$ with $r=0.999$ and $\sigma_{a}=9.8*10^{-3}$, i.e., as expected a slightly larger value than $2/7$. After breakup at $t_{+}$ we find $a=0.26$ with $0.9995$ and $\sigma_{a}=6.6*10^{-3}$. We remark that after the breakup time, the asymptotic solution is obtained by gluing the Hopf solution and the theta-functional solution (\[elliptic\]).

These results indicate that the scalings in (\[KdVrescaled\]) are indeed observed by the KdV solution. We show the corresponding situation for $t_{-}$ for two values of $\epsilon$ in Fig. \[figscalings\].

Multiscale solution
-------------------

In Fig. \[figpain4breakg\] we show the numerical solution of the KdV equation for the initial data $u_{0}$ and the corresponding PI2 solution (\[KdVrescaled\]) for $\epsilon=10^{-2}$ close to breakup. It can be seen that the PI2 solution (\[KdVrescaled\]) gives a correct description of the KdV solution close to the breakup point. For larger values of $|x-x_{c}|$ the multiscale solution is not a good approximation of the KdV solution.

A similar situation is shown in Fig. \[figpain4break1e6\] for the case $\epsilon=10^{-3}$. Obviously the approximation is better for smaller $\epsilon$. Notice that the asymptotic description is always better near the leading edge than near the trailing edge.

In Fig. \[figpain4delta\] we plot in green the difference between the PI2 multiscale solution and the KdV solution and in blue the difference between the KdV solution and the asymptotic solutions (\[Hopfsol\]) and (\[elliptic\]). It is thus possible to identify a zone around $x_{c}$ in which the multiscale solution gives a better asymptotic description. The limiting values of this zone rescaled by $x_{c}$ are shown in Fig. \[figwidthp12c\] for the critical time. It can be seen that the zone always extends much further to the left (the direction of propagation) than to the right.

The width of this zone scales roughly as $\epsilon^{3/7}$, more precisely we find ${\epsilon}^a$ with $a=0.468$, $r=0.981$ and $\sigma_{a}=0.073$. We observe that the numerical scaling is smaller than the one predicted by the formula (\[scalings\]). The matching of the multiscale and the Hopf solution can be seen in Fig. \[figdelta2ec\].

For larger times, the asymptotic solution (\[Hopfsol\]) and (\[elliptic\]) gives as expected a better description of the KdV solution, see Fig. \[figpain4delta1e6.226\] for $\epsilon=10^{-3}$ and $t=0.226$. Close to the leading edge, the oscillations are, however, better approximated by the multiscale solution.

To study the scaling of the difference between the KdV and the multiscale solution, we compute the $L_{\infty}$ norm of the difference between the solutions in the rescaled $x$-interval (\[interval\]) with $\alpha=3$. We find that this error scales at the critical time roughly like $\epsilon^{5/7}$. More precisely we find a scaling ${\epsilon}^a$ where $a=0.708$ ($5/7=0.7143\ldots$) with correlation coefficient $r=0.9998$ and standard error $\sigma_{a}=0.012$. Before breakup at the times $t_{-}(\epsilon)$ we obtain $a=0.748$ with $r=0.9996$ and $\sigma_{a}=0.016$, after breakup at the times $t_{+}(\epsilon)$ we get $a=0.712$ with $r=0.9999$ and $\sigma_{a}=6.2*10^{-3}$. Notice that the values for the scaling parameters are roughly independent of the precise value of the constant $\alpha$ which defines the length of the interval (\[interval\]). For instance for $\alpha=2$, we find within the observed accuracy the same value. In [@CV] Claeys and Vanlessen showed that the corrections to the multiscale solution appear in order $\epsilon^{3/7}$. For the values of $\epsilon$ we could study for our KdV example, the corrections are apparently of order $\epsilon^{5/7}$.

Outlook
=======

The Camassa-Holm equation [@CH] (see also [@Fo]) $$\label{CH}
u_t+6uu_x-\epsilon^2(u_{xxt}+4u_xu_{xx}+2uu_{xxx})=0$$ admits a bi-Hamiltonian description after the following Miura-type transformation $$\label{cam-holm1}
m=u-\epsilon^2 u_{xx}.$$ One of the Hamiltonian structure takes the form $$\label{cam-holm-pb1}
\{ m(x), m(y)\} =\delta'(x-y) -\epsilon^2 \delta'''(x-y)$$ so that the Camassa-Holm flow can be written in the form $$\label{hamCH}
m_{t}=\{m(x), H\},\quad H= \int (u^3+uu_x^2)dx.$$ To compare the Hamiltonian flow in (\[riem2\]) with the one given in (\[hamCH\]) one must first reduce the Poisson bracket to the standard form $\{ \tilde u(x), \tilde u( y)\}_1=\delta'(x-y)$ by the transformation $$\tilde u =\left( 1-\epsilon^2 \partial_x^2\right)^{-1/2} m=m+\frac12 \epsilon^2 m_{xx} +\frac38 \epsilon^4 m_{xxxx}+\dots.$$ After this transformation, the Camassa-Holm equation will take for terms up to order $\epsilon^{4}$ the form $$\tilde u_t +6\tilde u\, \tilde u_x +\epsilon^2 ( 8\tilde u_{x} \tilde u_{xx} + 4\tilde u\, \tilde u_{xxx}) +\epsilon^4( 20\, \tilde u_{xx} \tilde u_{xxx} + 12\, \tilde u_x \tilde u_{xxxx} +4\tilde u\, \tilde u_{xxxxx})+\dots=0 .$$ which is equivalent to (\[riem2\]) after the substitution $$c=48\tilde{u}, \quad p =2\tilde{u}.$$ At the critical point $(x_c,t_c,u_c)$ the Camassa-Holm solution behaves according to the conjecture in [@dubcr] as $$u(x,t,\epsilon) = u_c -\left(\dfrac{{\epsilon}^2|c_0|}{  k^2}\right)^{1/7} U(X,T)+O({\epsilon}^{\frac{4}{7}}),\quad c_0=4u_c$$ where $$X =-\dfrac{1}{\epsilon}\left(\dfrac{{\epsilon}}{k|c_0^3|}\right)^{1/7} 
( x - x_c - 6 u_c (t-t_c)),\quad
T=\left(\dfrac{1}{{\epsilon}^4k^3c_0^2}\right)^{1/7}(t-t_c)$$

In Fig. \[figpainch9\] we show the numerical solution to the CH equation for the initial data $u_{0}=-\mbox{sech}^{2}(x)$ and $\epsilon=10^{-2}$ at several values of time near the point of gradient catastrophe of the Hopf equation. It is interesting to compare this to the corresponding situation for the KdV equation in Fig. \[figpain4breakg\]. It can be seen that there are no oscillations of the CH equation on left side (the direction of the propagation) of the critical point, whereas in the KdV case all oscillations are on this side. The quality of the approximation of the CH and the KdV solution by the multiscale solution is also different. In the KdV case, the solution is well described by the multiscale solution on the leading part which includes the oscillations, whereas the approximation is less satisfactory on the trailing side. A similar behavior is observed in the CH case, but since the oscillations are now on the trailing side, they are not as well approximated as in the KdV case. The leading part of the solution near the critical point is, however, described in a better way.

The same qualitative behavior can also be observed for smaller $\epsilon$ in Fig. \[figpainch4\], though the quality of the approximation increases as expected on the respective scales. Note that we plotted in Fig. \[figpainch9\] and Fig. \[figpainch4\] the CH solution instead of the function $\tilde{u}$, since there are no visible differences between the two for the used values of $\epsilon$.

Numerical solution of the fourth order equation
===============================================

We are interested in the numerical solution of the fourth order ordinary equation (ODE) (\[PI2\]) with the asymptotic conditions (\[PI2asym\]). Numerically we will consider the equation on the finite interval $[X_{l},X_{r}]$, typically $X_{r}=-X_{l}=100$. In the exterior of this interval the solution to the equation (\[PI2\]) is obtained in the form of a Laurent expansion of $F$ around infinity in terms of $Y=X^{1/3}$, $$U = Y + \sum_{n=1}^{\infty}\frac{(-1)^{n}a_{n}}{Y^{n}}.
    \label{pain3}$$ We find the non-zero coefficients (not-given coefficients vanish) $a_{1}= 2T,$ $a_{5}= - 8T^3/3$, $a_{6}= 1/18$, $a_{7}=16T^4/3$, $a_{8}=- 5T/27$, $a_{10}=14T^2/27$, $a_{11}=- 256T^6/9$, $a_{12}=16T^3/3$, $a_{13}=640T^7/9-7/108$, …This expansion also determines the boundary values we impose at $X_{l}$, $X_{r}$ for $U$ and $U_{X}$.

The solution in the interval $[X_{l},X_{r}]$ is numerically obtained with a finite difference code based on a collocation method. The code *bvp4c* distributed with Matlab, see [@bvp4c] for details, uses cubic polynomials in between the collocation points. The ODE (\[PI2\]) is rewritten in the form of a first order system. With some initial guess (we use $U_{0}=-X^{1/3}$ as the initial guess), the differential equation is solved iteratively by linearization. The collocation points (we use up to 10000) are dynamically adjusted during the iteration. The iteration is stopped when the equation is satisfied at the collocation points with a prescribed relative accuracy, typically $10^{-6}$. The values of $U$ in between the collocation points are obtained via the cubic polynomials in terms of which the solution has been constructed. This interpolation leads to a loss in accuracy of roughly one order of magnitude with respect to the precision at the collocation points. To test this we determine the numerical solution via *bvp4c* for (\[PI2\]) on Chebychev collocation points and check the accuracy with which (\[PI2\]) is satisfied via Chebychev differentiation, see e.g. [@trefethen1]. We are interested here in values of $|T|<1$ and $|X|<10$. It is found that the numerical solution with a relative tolerance of $10^{-6}$ on the collocation points satisfies the ODE to the order of better than $10^{-4}$, see Fig. \[figpain4test\] where we show the residual $\Delta$ by plugging the numerical solution into the differential equation. It is straight forward to obtain higher accuracy by requiring a lower value for the relative tolerance, but we will only need an accuracy of the solution of the order of $10^{-4}$ here.

[99]{}

A. Bressan, *One dimensional hyperbolic systems of conservation laws*, Current developments in mathematics (2002), 1–37, Int. Press, Somerville, MA, 2003.

E. Brézin, E. Marinari and G. Parisi, *A nonperturbative ambiguity free solution of a string model*, Phys. Lett. B, [**2**42]{} (1990), no. 1, 35–38.

R. Camassa and D. D. Holm, An integrable shallow water equation with peaked solitons, Phys. Rev. Lett. [**71**]{} (1993), 1661-1664.

T. Claeys and M. Vanlessen, *The existence of a real pole-free solution of the fourth order analogue of the Painleve I equation*, Preprint:http://xxx.lanl.gov/math-ph/0604046.

T. Claeys and M. Vanlessen, *Universality of a double scaling limit near singular edge points in random matrix models*, Preprint:http://xxx.lanl.gov/math-ph/0607043. P. Deift, and X.  Zhou, *A steepest descent method for oscillatory Riemann-Hilbert problems. Asymptotics for the MKdV equation*, Ann. of Math. (2), **137**, (1993), 295–368.

P. Deift, S. Venakides, and X. Zhou, *New result in small dispersion KdV by an extension of the steepest descent method for Riemann-Hilbert problems*, IMRN [**6**]{}, (1997), 285-299.

P. Deift, T. Kriecherbauer, K. T.-R. McLaughlin, S. Venakides, and X. Zhou, *Uniform asymptotics for polynomials orthogonal with respect to varying exponential weights and applications to universality questions in random matrix theory*, Comm. Pure Appl. Math. [**5**2]{} (1999), no. 11, 1335–1425.

B. Dubrovin and S. P. Novikov, *A periodic problem for the Korteweg-de Vries and Sturm-Liouville equations. Their connection with algebraic geometry*. Dokl. Akad. Nauk SSSR , (1974), 531–534.

B. Dubrovin, *On Hamiltonian Perturbations of Hyperbolic Systems of Conservation Laws, II: Universality of Critical Behaviour*, Comm. Math. Phys., **267** (2006), 117.

B. Dubrovin and Y. Zhang, *Normal forms of hierarchies of integrable PDEs, Frobenius manifolds and Gromov - Witten invariants,* Preprint:http://xxx.lanl.gov/math.DG/0108160.

A. S. Fokas, *On a class of physically important integrable equations*, Physica [**D 87**]{} (1995), 145–150.

T. Grava and C. Klein, *Numerical solution of the small dispersion limit of Korteweg de Vries and Whitham equations*, to appear in Comm. Pure Appl. Math. (2006).

A. G. Gurevich and L. P. Pitaevskii, *Non stationary structure of a collisionless shock waves*, JEPT Letters [**1**7]{} (1973), 193-195.

A. Its and V. B. Matveev, *Hill operators with a finite number of lacunae*, (Russian) , Funkcional. Anal. i Priložen. [**7**]{} (1975), no. 1, 69–70.

S. Kamvissis, K.D.T.-R  McLaughlin, P. Miller, *Semiclassical soliton ensembles for the focusing nonlinear Schrödinger equation*, Annals of Mathematics Studies, **154**, Princeton University Press, Princeton, NJ, 2003.

C. Klein, *Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger equation*, preprint (2006).

P. D. Lax and C. D. Levermore, *The small dispersion limit of the Korteweg de Vries equation, I,II,III*, Comm. Pure Appl. Math. [**36**]{} (1983), 253-290, 571-593, 809-830.

L. F. Shampine, M. W. Reichelt and J. Kierzenka, *Solving Boundary Value Problems for Ordinary Differential Equations in MATLAB with bvp4c*, available at

        http://www.mathworks.com/bvp_tutorial
        

Fei-Ran Tian, *Oscillations of the zero dispersion limit of the Korteweg de Vries equations*, Comm. Pure App. Math. [**4**6]{} (1993) 1093-1129.

A. Tovbis, S. Venakides, X. Zhou,*On semiclassical (zero dispersion limit) solutions of the focusing nonlinear Schrödinger equation*. Comm. Pure Appl. Math. [**57**]{} (2004), no. 7, 877–985.

L. N. Trefethen, [ *S*pectral Methods in MATLAB]{}, SIAM, Philadelphia, PA, 2000.

S. Venakides, *The Korteweg de Vries equations with small dispersion: higher order Lax-Levermore theory*, Comm. Pure Appl. Math. [**43**]{} (1990), 335-361.

G. B. Whitham, [*L*inear and nonlinear waves]{}, J.Wiley, New York, 1974.

[^1]: We thank G. Carlet, B. Dubrovin and J. Frauendiener for helpful discussions and hints. CK and TG acknowledge support by the MISGAM program of the European Science Foundation. TG acknowledges support by the RTN ENIGMA and Italian COFIN 2004 “Geometric methods in the theory of nonlinear waves and their applications”.
