Peakons of Novikov Equation via the Homotopy Analysis Method

: The problem of the peakon and antipeakon solutions of the Novikov equation including the term ω u x is studied. It is well known that, when ω = 0, the Novikov equation admits peakon and antipeakon solutions. In this study, it is shown via the homotopy analysis method that, even in the case where ω (cid:54) = 0, the Novikov equation also admits peakon and antipeakon solutions.


Introduction
Both the celebrated KdV equation and the modified KdV (mKdV) equation where α, β, and γ constants, are among the simplest equations describing wave propagation. Equation (1) was originally introduced for describing propagation of unidirectional shallow water waves, whereas Equation (2) describes nonlinear acoustic waves. However, both of the aforementioned equations were also successfully used in several other physical problems; see, for example, [1] and the references therein. Several generalizations of them have appeared in the literature. One of them is the well known Camassa-Holm (CH) equation [2] u t + 2ωu x − u xxt + 3uu x = 2u x u xx + uu xxx , u = u(x, t) where ω is constant, which is also an equation describing dispersive shallow water wave propagation. Associated with the CH equation are also the Degasperis-Procesi (DP) equation and the Novikov equation (NE) u t − u xxt + 2ωu x = u 2 u xxx + 3uu x u xx − 4u 2 u x , u = u(x, t).
The DP equation for ω = 0 was derived in [3] in the context of a kind of asymptotic integrability of dispersive nonlinear equations, whereas the NE for ω = 0 was derived in [4] in an attempt to classify generalized CH type equations possessing infinite hierarchies of higher symmetry. Both Equations (4) and (5) can be thought of as generalizations of the CH equation. All Equations (1)-(5) share a common feature: they are nonlinear integrable equations with soliton solutions. However, among them, only the mKdV and the Novikov equations contain a cubic nonlinearity.
There is an extended literature for all Equations (1)- (5) where several mathematical issues for the behavior of their solutions are addressed. Apart from the pure mathematical interest on these equations, an important reason for studying them is the fact that they constitute simplified models for shallow water waves, which are governed by Euler equations, a set of equations quite difficult to be solved in general. Thus, even the study of the "simpler", compared to Euler, Equations (1)-(5), can provide useful information on the problem of formation and propagation of shallow water waves.
Equations (3)-(5) share one more common feature: they all have peakon solutions, i.e., peaked-shaped soliton solutions, which are solitons with discontinuous first derivative. Peakon solutions are weak solutions of (3)-(5), but they are orbitally stable. There are many reasons why peakons are interesting and some of them are discussed in [5,6], including the fact that they "are reminiscent of Stokes waves" for which "there is no closed form available and the peakons capture the essential features of these extreme waves." In addition, peakons are exact solutions that are always important when studying differential equations. Peakons of the CH equation were already presented in [2] and are studied for all Equations (3)-(5), mostly from a theoretical point of view. This paper will be focused on the Novikov Equation (5), for which there exist several papers where various theoretical aspects for (5) are studied. Indicatively, a few recent works for the case ω = 0 will be briefly mentioned here. In [7], a matrix Lax pair as well as a bi-Hamiltonian structure for NE were found. In [8], modern group analysis was utilized in order to find the Lie point symmetries of NE and to prove that it is self-adjoint. The stability of peakons of NE was studied in [6,9]. The periodic and non-periodic Cauchy problems for NE were studied in [10][11][12][13][14][15] with respect to its well-posedness, dependence on initial data, continuity results for the data-to-solution map, existence, and uniqueness of global strong and weak solutions in Sobolev and Besov spaces, existence, and uniqueness of real analytic solutions, as well as persistence properties of the strong solution. In [16], sufficient conditions on the initial data were established in order to guarantee the formation of singularities in finite time of the solutions of NE. In [5], explicit multipeakon solutions of NE were calculated by inverse spectral methods using the matrix Lax pair found in [7]. In [17], explicit formulas were given for the characteristic curves of the CH, DP, and Novikov equations. In particular, Section 2 of [17] consists of a very good survey on the existing knowledge on these three equations. In [18], it was proved that a generalized Camassa-Holm equation, which includes the CH, the DP, and the Novikov equations, possesses peakon traveling wave solutions in Sobolev spaces. In [19], 2-peakon solutions on both the line and the circle for NE were constructed, with an asymmetric antipeakonpeakon initial profile. The peakon and antipeakon both move in the positive direction and collide at a finite time resulting in ill-posedness of the corresponding Cauchy problem. The introduction of [19] summarizes well the existing knowledge on the properties of NE.
Since NE is a nonlinear equation, it is not easy to explicitly find exact solutions of it, especially peakons, via analytic methods. Thus, numerical or approximate methods are usually implemented. Although there exists an excessively large amount of papers which numerically investigate smooth solitons of various well-known partial differential equations (PDEs), there are much less papers where peakons are numerically computed. Even more, in these papers, more sophisticated techniques than the ones used for numerically computing smooth solitons are implemented such as a finite volume method [20], a pseudospectral scheme [21], or a local discontinuous Galerkin method [22].
Another choice for finding and studying exact solutions of differential equations, including solitons and peakons, are via approximate techniques with the aid of which it is possible to analytically calculate an approximation of the corresponding solution. Nowadays, the rapid development and use of symbolic computation software has greatly aided the implementation of approximate methods. Among the most popular methods are perturbation methods (see, for example, [23] or [24]), the Adomian decomposition method (ADM) [25], the differential transform method (DTM) [26], and the homotopy analysis method (HAM), which is described in detail and applied to a variety of physical and mathematical problems in the books [27,28]. It is worth mentioning that the DTM, combined with the method of steps, has been successfully applied for solving PDEs [29] and boundary value problems for delayed differential equations [30].
Among the aforementioned approximate techniques, one of the most popular, and, at the same time robust method, is HAM, which, since its appearance, has been extensively used in a variety of physical problems, such as beam problems, oscillation problems, fluid mechanics problems, and, of course, wave propagation problems. Probably the first papers where HAM was applied to the study of nonlinear waves were [31,32], in which solutions to the problem of progressive waves in deep water and solutions for periodic waves for the mKDV equation, respectively, were derived. Ever since, HAM has been successfully applied to obtain soliton and/or peakon solutions of various important nonlinear PDEs including but not limited to the Fitzhugh-Nagumo equation [33], the CH equation [34,35], a 5th order KdV equation [36], and the DP equation [37,38].
Applying HAM to problems of nonlinear water waves, including the CH equation, in order to obtain peakon solutions, is one of the five challenging problems included in [39]. Motivated by this and the importance of NE and its peakons, the aim of the present paper is to apply HAM in order to obtain peakon and antipeakon solutions of NE. The mathematical formulation of the problem, as well as a sort description of HAM applied to NE, are given in Section 2. The implementation of HAM was achieved using Mathematica. In Section 3, a thorough discussion of the obtained results is given with respect to several values of ω, the velocity of the peakon, and initial guesses of the solutions to be obtained. Moreover, the corresponding averaged squared residual errors are calculated. The results for NE are similar to the results obtained for DP by Liao in [38], which also inspired this work. Finally, in Section 4, some brief conclusions are deduced.

Mathematical Formulation and Implementation of HAM
It is well known that traveling wave solutions of a PDE are obtained via the transformation where c = 0 is the velocity of the traveling wave. Solitons in particular of a PDE are solutions of the form (6) satisfying the boundary conditions g, g , g , . . . → 0, when η → ±∞.
Thus, solitons of NE are obtained as solutions of the following boundary value problem (BVP) g, g , g → 0, when η → ±∞, which is obtained after substituting (6) into (5) and taking into consideration (7). Of course, denotes differentiation with respect to η.
Although η ∈ R, it suffices to study (8) and (9) only for η ≥ 0, since the simple change of variables y = −η, g(η) = f (y), leaves (8) unchanged. Moreover, it is easy to see that, if g is a solution of (8) and (9), so is −g. Thus, if a peakon is a solution of NE, so is its corresponding antipeakon.
It is obvious that, when ω = 0, the wave speed c does not appear in (13).
In order to study (13) and (14), the HAM will be applied. According to HAM, the corresponding solution is given in the form of the so-called homotopy series where w 0 (η) is an initial approximation of the solution of (13) and (14). The functions w m (η) are the solutions of the linear BVPs where L is a linear differential operator associated with (13),h is a nonzero auxiliary quantity called convergence control parameter, H(η) is a nonzero auxiliary function and with w s = w s (η), s = 1, . . . , m − 1. All of w 0 (η), L, and H(η) are appropriately chosen in such a way so that w m (η), m = 1, 2, . . . and consequently the series (15) are expressed only in terms of a proper set of base functions adequately describing the solution of the problem under consideration. Due to the form of (11) and the fact that exponential functions have already been successfully used for the derivation of peakons via HAM, an appropriate set of base functions for obtaining peakons of NE is the set where κ = 1 − 2ω c . Obviously, for κ to be real, the inequality should hold. It should also be mentioned that the corresponding to (13) linear equation has exponential solutions only in the case where (21) holds.
Thus, the solution of (13) and (14) will be expressed in the form where the coefficients a m are obtained from (23) after substituting w m (η). In view of the above, an appropriate choice for w 0 (η), L and H(η) is the following: where A a nonzero parameter. In practice, M terms of the homotopy series are calculated, in order to obtain an adequate approximation of w(η). The corresponding series is called an approximate solution of M−order. The convergence control parameterh will be chosen in such a way so that the series on the right-hand side of (15) is convergent. When this is true, it is known (see [27] or [28]) that the solution found by HAM is a solution of (13) and (14). One way to chooseh is to use the so-calledh−curves, which were proposed in [27] and have been greatly used. Theh−curves are graphs of specific, usually important, quantities associated with the solution of the DE under consideration, calculated at specific values of the independent variable, againsth. If (15) is convergent, a horizontal line segment will be apparent in theh−curves. Furthermore, if a value ofh is chosen from this region, it is quite certain that (15) converges.
Another way to chooseh is to find the optimumh for which the averaged squared residual error is defined by where is minimized, see [28] [Chapter 3]. In the present paper, the second approach will be used for reasons explained in Section 3.

Results and Discussion
In the current section, the solution of (13) and (14) will be determined using the homotopy analysis method, for various values of ω, c, and A. The implementation of HAM has been carried out via Mathematica, following the outline given in Section 2 for the nonlinear operator R m defined by (19) and the choices of w 0 (η), H(η) and L given in (24).
To begin with, the case where c = 1, ω = 0.1, and A = 0.1 is considered. In order to calculate the solution and at the same time find an appropriate value ofh, the 10th order approximation is calculated. Due to the geometry of the solution, a quantity of interest for w(η) is w (0 + ) = lim η→0 + w (η). Thus, the correspondingh-curve is graphed in Figure 1. It is obvious from this graph that a valid region forh is, roughly speaking, from −0.5 up to −2. Thus, it is natural to chooseh = −1.5. Table 1 shows w (0 + ) as well as the corresponding averaged squared residual error defined by (26), forh = −1.5 and for several values of the order M from 2 up to 16. For the calculation of the errors E M , it suffices to choose a = 10, since the solution of (13) and (14) decay exponentially to zero.  Table 1 that the errors E M decay quickly to order 10 −19 and the value of w (0 + ) converges to −0.0893676. However, as mentioned in [28], the whole procedure can be highly facilitated and accelerated by finding the optimumh of a specific approximation, via minimizing the corresponding errors and then use this optimumh in order to calculate higher approximations of the solution. Thus, in the rest of the paper, this approach will be used. The optimumh for the 6th order approximation is found to be −1.00791 and, by using the rounded valueh = −1, the 10th order approximation of the solution is calculated. This peakon solution, together with its corresponding antipeakon solution, is presented in Figure 2. The second row of Table 2 presents the values of w (0 + ) and E M in the case where c = 1, ω = 0.1, A = 0.1, andh = −1. It is obvious that the value −0.0893676 of w (0 + ) is attained at the second approximation and that the errors are much smaller than in the case whereh = −1.5, reaching the order 10 −47 at the 10th approximation.

It is obvious from
At the same table, Table 2, the values of w (0 + ) and E M are also presented for greater values of A, namely for A = 1/5, A = 3/10, and A = 2/5. It is obvious that the HAM converges rapidly again and the corresponding errors are very small. However, the obtained results indicate that, as A increases, the errors for the same 10th approximation of the solution of (13) and (14) gradually increase, although they still remain very small. The corresponding peakon solutions for c = 1, ω = 0.1,h = −1 and the aforementioned values of A are shown in Figure 3.  Table 2. w (0 + ) and corresponding errors for c = 1, ω = 0.1,h = −1 and A = 1/10, 1/5, 3/10, and 2/5. Similarly, in Table 3, the values of w (0 + ) and E M are presented again forh = −1, in the case where c = 1 and A = 1/10, but for four different values of ω, namely for ω = 0.1, ω = 0.2, ω = 0.3, and ω = 0.4. It is obvious again that HAM rapidly converges, but in contrast with the results of Table 2, the errors don't seem to be affected by increasing ω. These peakon solutions are presented in Figure 4. The results indicate that the solution decays more slowly, as ω increases.
Several simulations have been carried out also for larger values of c and ω. The corresponding results indicate similar conclusions as before. The peakon solutions of problem (13) and (14) for c = 4 and c = 36 and for various values of ω are presented in Figures 5 and 6, respectively. Table 3. w (0 + ) and corresponding errors for c = 1, A = 1/10,h = −1 and ω = 0.1, 0.2, 0.3 and 0.4.

Conclusions
It is well known that the Novikov Equation (5) has peakon solutions in the case ω = 0. In this paper, it was shown via the homotopy analysis method that (5) also possesses peakon solutions even in the case where ω = 0. In order to apply HAM, an initial approximation of the solution must be used, which in this case was the function w 0 (η) = Ae −κ|η| , with κ = 1 − 2ω c . Several analytic approximations of the peakon and antipeakon solutions of (5) were obtained via HAM for various values of ω, A and the velocity c. The results show that the solutions decay more slowly, as ω increases. Moreover, the averaged squared residual errors were calculated for the aforementioned approximations. The results indicate that these errors are not affected when ω increases, but are increased when A increases. The obtained results in this paper for NE are similar to the results of [38] for the DP equation.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: