Camassa and Holm  [4] introduce the equation,

discuss its analytical properties and sketch its derivation. Here $u$ is the fluid velocity in the $x$ direction, $\kappa$ is a constant related to the critical shallow-water wave speed, and subscripts denote partial derivatives.

Actually, the equation is a new completely integrable dispersive shallow water equation that is bi-Hamiltonian. For $\kappa :=0$ the authors quoted above showed that this equation has solitary water waves of the form $c\exp \left(-|x-ct|\right)$ which they called ”peakons”. Dropping both high-order terms from the r. h. s. of (1) leads to the familiar Benjamin-Bona-Mahony (BBM) equation  [2], or at the same order, to the Korteweg-de Vries (KdV) equation. Nevertheless, these extra terms profoundly alter the behavior of solitary waves. Thus, this extension of the BBM equation possesses soliton solutions whose limiting forms as $\kappa \to 0$ have peaks where first derivatives are discontinuous, i.e. peakons. The lack of smoothness at the peak of the peakon introduces high-frequency dispersive errors into the calculation. This is one of the main numerical challenges.

There is a fairly rich literature gathered around the Camassa-Holm equation. For instance, in  [3] J. P. Boyd derives a perturbation series solution for general $\kappa$ which converges even at peakon limit and also gives three analytical formulations for spatially representations of the peakons. He also observes that although the Camassa-Holm equation is integrable for general $\kappa$, it appears that imbricating solitary waves generates an exact spatially periodic solution only for the special cases $\kappa =0$ and $\kappa /c=1/2$. Camassa, Holm and Hyman [5] present periodic numerical solutions to equation (1) and a large discussion to this equation as a Hamiltonian system. They focus on the case $\kappa =0$. Actually, we also have some performing methods, like Runge-Kutta discontinuous Galerkin, for periodic case or bounded domains obtained by empiric truncation.

The main motivation to design this MoL-SiC scheme is to avoid the arbitrary truncation of the domain or periodicity assumptions. Instead, we suggest an asymptotic truncation of the integration domain through the scaling factor of SiC. The present paper provides accurate numerical results for arbitrary $\kappa$ case. We use a MoL based on SiC in order to discretize the spatial derivatives and TR-BDF2 in order to march in time.

The boundary conditions for solution $u(x,t)$ will be taken to be regular at infinity, i.e.

for at least $n=0,1,2$ where $u_x^{\left(n\right)}\left(x\right)$ denotes the $nth$ derivative with respect to $x$. They are here behavioral (natural) rather than numerical and have been inspired by the analysis in [2].

The pure initial-value problem for the above equation is to ask for a solution $u$ defined for $(x,t)\in ℝ\times {ℝ}^{+}$ having a specified initial configuration namely

A typical situation arises wherein the initial disturbance has sensibly finite extent. Thus, we suppose a datum $u_0\left(x\right)$ being drawn from some classes of functions having limit at $\pm \mathrm{\infty }$, i.e., satisfying at least the first condition in (2). In such circumstances the pure initial value problem for (1) is well-posed in Hadamard’s classical sense (see [1]). That is, corresponding to a suitable initial datum $u_0\left(x\right)$, there is a unique function $u(x,t)$ defined on $ℝ\times {ℝ}^{+}$ satisfying the differential equation there, and for which (2) holds. Moreover, if $u_0\left(x\right)$ is slightly perturbed in its function class, then the solution $u$ changes only slightly in response.

The MoL solution to the equation (1) reads

where $S(j,h)\left(x\right)$ are the sinc discrete functions introduced and analyzed by F. Stenger in [10] and $u_j\left(t\right)$ for $j:=1,\mathrm{\dots },N$ are the nodal time unknowns. $N\in \mathbb{N}$ is the order of approximation.

The functions $S(j,h)\left(x\right)$ are defined by

where $x_j$ are the equispaced nodes with spacing $h$, symmetric with respect to the origin and ${𝐱}_{Si}:={(x_1,\mathrm{\dots },x_N)}^T$.

It is important to observe that SiC approximation $u_N(x,t)$ defined by (4) satisfies the regularity condition at infinity (2). The sinc differentiation collocation matrices on ${𝐱}_{Si}$ are denoted by ${𝐃}_{Si}^{\left(n\right)}$, $n=1,2,3$. In order to implement the MoL we have used the form of these matrices provided by Weideman and Reddy in  [12].

The Camassa-Holm PDE is one of the form

where the operators $ℳ$ and $ℒ$ are linear and $𝒩$ is a nonlinear one. Once we discretize the spatial part of this PDE we get a system of ODEs,

where the vector $𝐔$ is defined by $𝐔:={(u_1\left(t\right),\mathrm{\dots },u_N\left(t\right))}^T.$

The particular form of the discrete formulation (6) for Camassa-Holm equation reads

where the dot, as in MATLAB, signifies the element wise product of two vectors.

For our dispersive problem it is known that the eigenvalues of the linear operator are close to the imaginary axis (see the analysis in [9]). They are depicted in Fig. 2 and are situated inside or very close to the region of stability of TR-BDF2 depicted in Fig. 1. This is the main reason in choosing this finite difference scheme to solve the ODEs system (7).

We adopt the strategy from [11] in order to establish the linear (numerical) stability of MoL. We have successfully used this strategy in our contribution [7] in solving numerically the BBM equation which is another pseudo-parabolic one. The linearization of the discrete equation (7) is simply

where $𝐈$ is the unit matrix of order $N$. Its pseudospectrum is depicted in Fig. 2 and shows that indeed all eigenvalues are on the imaginary axis.

If we overlap Fig. 2 and the region of stability of TR-BDF2 from Fig. 1 we observe the $\epsilon$-pseudospectrum lies within a distance $O\left(\epsilon \right)$+$O\left(\mathrm{\Delta }t\right)$ of the stability region as $\epsilon \to 0$ and $\mathrm{\Delta }t\to 0$. Typically the time step has been of order ${10}^{-2}$ when (1) has been solved.

Thus we can conclude that the MoL-SiC scheme along with TR-BDF2 is numerically stable irrespective of values of parameter $\kappa \ne 0$. In fact, our numerical experiments have provided reasonable numerical results for the case $\kappa =0$ even if the above analysis of stability does not literally apply to this case.

Some numerical experiments will be reported in this section. We consider the interaction of two peakons which emerge from the initial data

Actually we solve the Cauchy’s problem (1)-(9). The interaction of two peakons as singular soliton solutions is depicted in Fig. 3.

The evolution of initial data (9) to its final form of two separated and reversed in order peakons is provided in Fig. 4.

These solitons travel with a speed proportional to their height and remain coherent after their collision (a short animation is instructive in this respect). In spite of the fact that peakons are nonanalytic solitons, a careful inspection of Fig. 3 and Fig. 4 shows that no oscillations (like Gibbs’ ones) appear when resolving the cusp of peakons.

The MoL-SiC along with TR-BDF2 proved to be a reliable and effective method to solve some Cauchy’s problems attached to a pseudo-parabolic equation. The computational effort is kept in fairly reasonable limits, i.e., some seconds in order to solve the system of semi-discrete ODEs on a usual machine.

The linear stability of MoL is established but a nonlinear one is an open problem. The conservation of two invariants is kept in reasonable limits during the time integration interval. The first one behaves better than the second as expected. Moreover, with respect to the accuracy of method we establish that the convergence is only exponentially.