Constructions of the soliton solutions to the good Boussinesq equation

The principal objective of the present paper is to manifest the exact traveling wave and numerical solutions of the good Boussinesq (GB) equation by employing He’s semiinverse process and moving mesh approaches. We present the achieved exact results in the form of hyperbolic trigonometric functions. We test the stability of the exact results. We discretize the GB equation using the finite-difference method. We also investigate the accuracy and stability of the used numerical scheme. We sketch some 2D and 3D surfaces for some recorded results. We theoretically and graphically report numerical comparisons with exact traveling wave solutions. We measure the L2 error to show the accuracy of the used numerical technique. We can conclude that the novel techniques deliver improved solution stability and accuracy. They are reliable and effective in extracting some new soliton solutions for some nonlinear partial differential equations (NLPDEs).


Introduction
The soliton theory is an important tool in describing various phenomena of wave propagation. Recently, several nonlinear equations have been successfully developed to investigate wave propagation. For instance, the Korteweg-de Vries (KdV) equation is mainly used to investigate the propagation of shallow water waves in one dimension, whereas the bad Boussinesq (BSQ) equation describes the wave propagation in two dimensions. The BSQ equation is also used in studying different processes appearing in electromagnetic waves in dielectrics [1], magnetosound waves in plasma [2], and the magnetoelastic waves in antiferromagnetic [3]. These equations and others have attracted the attention of a massive number of mathematicians and physicists since the 1970s due to their use in revealing the internal mechanisms of some sophisticated natural phenomena. Therefore numerous scientists have invented and discovered a wide variety of new methods to develop the soliton solutions for most NLPDEs. Some developed techniques and principles include the inverse scattering transform [4], the trial function process [5], the sine-cosine principle [6], the Weierstrass elliptic function approach [7], the tanh-sech technique [8], the Fexpansion technique [9], Hirota's bilinear principle [10], the modified tanh-function tech-nique [11], the extended tanh-technique [12], the exp(-f (ζ ))-expansion process [13,14], and the truncated Painleve expansion [15]. More information about other techniques can be obtained in [16][17][18][19][20][21][22][23][24][25][26][27].
The Boussinesq equation is given by where u(x, t) is a real-valued function, and subscripts represent partial derivatives. This equation models the propagation of shallow water waves in both directions. If the sign associated to u xxxx is changed, then we end up with a linearly stable equation, and the numerical computation can be achieved. Therefore the good Boussinesq (GB) equation is given by According to [28], the solitary waves described by the GB equation solely occur for a finite range of velocities and can merge into one solitary wave. The Boussinesq equation has been exactly and numerically solved using different approaches. For example, the modified decomposition method is applied in [29] to develop soliton solutions and periodic solutions, whereas a simplified version of the Hirota technique is used in [30] to extract several exact solutions for the GB problem. Furthermore, in [31], N-soliton solutions are determined by using the bilinear form. Nguyen [32] has employed the Hirota direct bilinear process to find the soliton solution of the GB problem. The power-law nonlinearity [33] and the variational iteration techniques [34] are derived to obtain solitary wave solutions for the Boussinesq equation and the GB equation, respectively. On the other hand, the GB equation is solved numerically by Ismail and Bratsos [35]. They have presented a conditionally stable scheme that is second-order in space and fourth-order in time. The authors in [36] have applied a finite difference process to the GB equation and discovered that the used technique is unconditionally stable and is fourth-order in space and second-order in time. Other numerical methods can be found in [8,[37][38][39][40][41][42][43][44][45].
In this work, we employ the He semiinverse approach and the moving mesh process [46] to construct the soliton and numerical solutions, respectively, for the GB equation. The numerical method followed in this paper uses a monitor function and moving mesh partial differential equations (MMPDEs) to distribute the grid of the points in the area where the error is high. The main advantage of this approach is the reduction of the error, especially in the curvature and variation regions, by distributing more points in such regions.
The boundary constraints are graphically deduced from the behavior of the exact traveling wave solutions as the time changes. The exact solutions vanish at the boundaries of the physical domain. As a result, the relevant boundary conditions are constants. In other words, u x = 0 and u xxx = 0 as x → ±∞.

Traveling wave solution
In this part, we focus on extracting the exact solution of the GB equation using the He semiinverse process [47,48]. Plugging the traveling wave transformation were w plays the role of the speed wave, into Eq. (2) and integrating twice with respect to ξ give Hence from Eq. (3) we obtain the following variational formulation: Next, we utilize a Ritz-like process to invoke bright optical solitons in the following form: where β and α are constants to be found later. Plugging Eq. (5) into Eq. (4) leads to To find the stationary points (fixed points) β and α of J, we differentiate J with respect to β and α and equate the results to zero: The solutions of this system are as follows: Substituting these values of β and α into Eq. (5) yields Hence the exact solutions of Eq. (2) can be simply expressed as where w indicates the speed parameter, and x 0 is a constant.

Stability analysis
We study the stability of the accomplished exact traveling wave results by employing the Hamiltonian system expression given by where w plays the role of the speed parameter, ϕ(w) is the momentum, and ρ indicates the obtained exact solution of Eq. (3). A sufficient condition for testing the stability is given by Applying the Hamiltonian system to Eq. (6) gives Taking the derivative with respect to w yields Hence we clearly observe that the exact solution is stable for all ξ ∈ (-∞, +∞) under the parameter value |w| < 1. The sign of w plays a crucial role in determining the direction of the time evolution of the traveling wave solution. Consequently, we take w = 0.5 and x 0 = -10 to plot the following figures.

Numerical investigation
In this section, we numerically solve GB equation via the adaptive moving mesh method. This novel approach uses a monitor function, which distributes the mesh nodes along the evolving solution at each time step. To discretize in space, we utilize the central finite differences. To study Eq. (2) numerically, we assume that Differentiating Eq. (7) with respect to t and substituting the result into Eq. (2) lead to Hence Eq. (2) is converted to the following system: The relevant boundary conditions are obtained from Fig. 1, from which we observe that where a and b are the endpoints of the physical domain. To find the numerical solutions of system (8)  Hence the exact solution of v(x, t) is given by Thus the initial conditions are given by Applying the adaptive moving mesh process to build the numerical solution of the GB problem requires equal subintervals for the movement of the mesh. Thus the physical domain [a, b] is first transformed into a computational domain, which is [0, 1], using the following transformation: Employing the physical coordinate x and the computational coordinate ζ leads to Thus the moving mesh connecting with the solutions u and v is formed as and the computational domain is formed as Since the solution u and the mesh x are functions of the variable t, applying the chain rule to system (8), we have In the moving mesh method the error indicator is selected to redistribute the mesh nodes where the solution presents large curvatures or variations. We attempted different MM-PDEs and monitor functions and we found that the obtained results are similar. As a result, we use the MMPDE5 [46,49,50] given by where τ ∈ (0, 1) is called a relaxation parameter, and (u, v) denotes the monitor function. A successful selection for this function often gives effective and dependable results. This can be attributed to the fact that the monitor function feeds the bending and variation regions in the solution with more nodes. Therefore we apply an exceptional function given by The discretization of the mesh Eq. (12) is given by subject to the boundary conditions The initial condition is given by splitting the physical domain into N x identical intervals as follows: where x j = a + j x with x = b-a N x and j = 1, 2, . . . , N x -1. Moreover, the discretizations of the coupled equations (11) are given by subject to the boundary conditions It is to be highlighted that the initial conditions are chosen by Eqs. (10) at t = 0. Here t presents the step size of time t. To evaluate the boundaries of u xx and v xx , we require some fictitious points Here we list the procedure of the alternating solution as follows: 1. At the time step t n the monitor function Eq. (13) is computed using the solutions u n and v n . 2. A new mesh at the time step t n+1 is obtained by solving scheme (14) using the monitor function n . The monitor function is fixed during the computation. 3. The solutions u n+1 and v n+1 are obtained by solving the numerical schemes (15) and (16) simultaneously using x n and x n+1 . 4. This procedure is repeated from step 1.

Accuracy of the numerical schemes
This subsection is devoted to introducing the accuracy of the adaptive mesh schemes (14), (15) and scheme (16). The approximated solutions u n j and v n j are replaced throughout by u(x j , t n ) and v(x j , t n ), respectively, at the point (x j , t n ). Then we introduce the step size for both the time and spatial variables. They are presented by t = t n+1t n , x + = x j+1x j , and x -= x jx j-1 . Taylor series expansion is used to determine x n+1 j+1 and x n+1 j-1 as follows: Adding these expressions and simplifying, we get In the same manner, we have Consequently, the truncation error of the implicit scheme (14) is To investigate the accuracy of schemes (15) and (16), we use Taylor series expansions for u n+1 j and u n j with respect to t/2 as follows: Subtracting these expansions and simplifying lead to The spatial first and second derivatives for u are approximated by the average of finite differences at t n and t n+1 as follows: Similarly, we can construct the approximation of u x and u xx at the time level t n . In the same manner, we can derive the second derivative for v at t n and t n+1 . Hence the truncation error of the adaptive moving mesh schemes (15) and (16) is given by As a result, the accuracy of the adaptive mesh scheme is of the second order in time and roughly more than the second order in space when x + = x -, that is,

Stability of the numerical schemes
Here we use the von Neumann analysis to examine the stability of the numerical scheme (14). We ignore the boundary conditions and consider (t n , ζ j ) with t n = n t, ζ j = j ζ , and x j = j x. Rearranging Eq. (14) yields where γ = t τ ζ 2 . To employ the von Neumann technique, we assume that where k is a constant. Substituting Eq. (20) into Eq. (19) gives Since (e ik ζ -2 + e -ik ζ ) = -4 sin 2 (k ζ /2), we have Hence Eq. (14) is unconditionally stable for γ ≥ 0. Now, to study the stability of scheme (15) and scheme (16), we assume that the mesh is fixed. Then the schemes are expressed by where δ 2 x = u n j+1 -2u n j + u n j-1 , μ = t/ x 2 , and α = t max 0≤j≤N x (1 + u n j ). We assume that u n j = λ n e ik xj , v n j = β n e ik xj .

Solving this system gives
Since the required conditions for the stability of the numerical schemes are |λ| ≤ 1, |β| ≤ 1, sin 2 ( 1 2 k x) ≤ 1, and α = t max 0≤j≤N x (1 + u n j ), the schemes are unconditionally stable for μ ≥ 0. In this work, we used MATLAB software for running the numerical schemes to obtain numerical results.    Fig. 3(a). The achieved solutions are  also shown similarly in 2D figures given in Figs. 4(a), (b). This leads to the conclusion that the used methods are applicable, powerful, and useful in solving other nonlinear PDEs. Table 1 presents a brief description of L 2 errors and CPU times consumed to arrive at t = 5 for the adaptive technique. Various values of x were used. The error decreases to reach 4.715 × 10 -6 when x = 1.0 × 10 -2 for a long time (0.819 × 10 +1 s). However, the method gives a small and acceptable error 2.19 × 10 -5 in short time 7.7 × 10 -1 when x = 2.5 × 10 -2 . The increase in the CPU time is due to the additional functions, which altogether evaluated along with the PDE. We can observe from Fig. 5 that the L 2 error dramatically decreases when the number of mesh points increases. However, the CPU time slightly increases. As a result, the adaptive moving mesh method is reliable and more computationally effective.

Discussion and results
In this paper, we successfully employed a novel technique to construct the numerical solutions of Eq. (2). This method distributes the mesh points into the regions in which the error is high. Taylor expansion was specifically applied on Eq. (14). The accuracy is of the first order in time and of the second order in space, as shown in Eq. (17). Furthermore, the accuracy of schemes (15) and (16) is of the second order in time and of slightly more than the second order in space when x + = x -, as described in Eq. (18). Moreover, the von Neumann stability showed that for Eq. (14), schemed (15) and (16) are unconditionally stable.
It is well known that the numerical solutions of PDEs are established by approximating the relevant equation on a mesh. Some solutions may contain dramatic spatial change. Hence a fine number of nodes is required on a small enough domain. However, this procedure is often extensive in computation. Consequently, we employed the adaptive moving mesh method, which sends more points into the regions in which there is a variation. The authors in [36] used the finite difference method to extract the numerical solution of the GB equation. The points are distributed on a uniform mesh. However, the adaptive moving mesh method gave new and more general results than those presented in [36]. We obtained more successful results during an acceptable time.

Conclusions
In this work, we discussed the traveling wave solutions and the numerical results of the good Boussinesq equation using the He semiinverse and the adaptive moving mesh methods, respectively. The accomplished exact solution was determined in the form of a hyperbolic trigonometric function. The exact solution was stable for all x ∈ (-∞, +∞) under the parameter values -1 < w < 1. The exact and numerical solutions were compared in some 3D figures, which showed that the solutions were almost the same. The stability and accuracy of the numerical scheme were studied. The L 2 error was measured to illustrate the accuracy of the adaptive moving mesh approach. The error vanished for a large number of mesh nodes, as shown in Fig. 5. The used methods gave successful and effective results for most nonlinear partial differential equations.