Research Article  Open Access
Solitary Wave Formation from a Generalized Rosenau Equation
Abstract
A generalized viscous Rosenau equation containing linear and nonlinear advective terms and mixed third and fifthorder derivatives is studied numerically by means of an implicit secondorder accurate method in time that treats the first, second, and fourthorder spatial derivatives as unknown and discretizes them by means of threepoint, fourthorder accurate, compact finite differences. It is shown that the effect of the viscosity is to decrease the amplitude, curve the wave trajectory, and increase the number and width of the waves that emerge from an initial Gaussian condition, whereas the linear convective term pushes the wave front towards the downstream boundary. It is also shown that the effect of the nonlinear convective term is to increase the steepness of the leading wave front and the number of sawtooth waves that are generated behind it, while that of the first dispersive term is to increase the number of waves that break up from the initial condition as the coefficient that characterizes this term is decreased. It is also shown that, for reasons of stability, the second dispersion coefficient must be much smaller than the first one and its effects on wave propagation are relatively small.
1. Introduction
In his 1986 and 1988 papers, Rosenau [1, 2] developed a formalism to treat the dynamics of discrete dense systems that can deal with wavewave and wavewall interactions that cannot be treated with the Kortewegde Vries (KdV) equation. Such a formalism leads to a quasicontinuum which is endowed with the leading effects due to discreteness.
Rosenau’s original equation may be written as where the subscripts and denote partial differentiation with respect to the time and the spatial coordinate, respectively. In the absence of the fifthorder derivative term, (1) reduces to a firstorder nonlinear wave equation which has analytical solutions and may have shock wave solutions if .
The Rosenau equation has been the subject of several analytical and numerical studies. Barreto et al. [3] prove the existence of solutions of (1) with the plus sign in the advectionlike term in moving domains by making use of the Galerkin method, multiplier techniques, and energy estimates. They also proved analogous results for the BenjaminBonaMahony or regularized longwave (RLW) equation with a linear advective term; that is, which models long waves in a nonlinear dispersive medium.
A onedimensional generalization of the Rosenau equation may be written as where denotes a forcing function and
Equation (4) with and corresponds to the original Rosenau equation and includes the linear (firstorder) wave equation, the inviscid and viscous Burgers equations, the RLW, modified RLW, and generalized RLW equations, the Kortewegde Vries (KdV) equation, the RosenauHyman equation, the CamassaHolm equation, the OlverRosenau equation, the Kawahara equation, the CooperShepherdSodano equation, and combinations thereof; these equations are presented in Table 1.

Park [4] studied the global existence and uniqueness of solutions of the following multidimensional generalized Rosenau equation:where is the Laplacian or harmonic operator, while the same author [5] provided pointwise decay estimates for (3) with , that is, a RosenauBurgers equation. Wang and Xu [6] studied the existence and global solution of the secondorder in time Rosenau equation and proved the existence of finitetime blowup for for and , by means of a potential well technique, while H. Wang and S. Wang [7] studied the longtime behavior of small solutions of the Cauchy problem for a (secondorder in time) Rosenau equation, obtained its global small solution, and analyzed the decay and scattering of such a solution. A similar study for the (firstorder in time) RosenauBurgers (RB) equation was performed by Liu and Mei [8] who proved that the solution of a nonlinear parabolic equation is a better asymptotic profile of the RB equation. H. Wang and S. Wang [9] proved the global existence of the solution to the Cauchy problem for the th dimensional Rosenau equation with a damping term when the initial data are small, while Kim and Lee [10] analyzed the convergence of a semidiscretization of .
Esfahani [11, 12] obtained solitary wave solutions to the generalized RosenauKdV and RosenauRLW equations, respectively, by means of the sech and trigonometric function methods, respectively, while Razborova et al. [13] studied the solitons, shock waves, and conservation properties of RosenauKdVRLW equations with power nonlinearities by means of a semiinverse variational method. Razborova et al. [14] used a perturbation method to study shallow water waves governed by the same equation. Choo et al. [15] obtained a posteriori error estimates of the Rosenau equation by means of a discontinuous Galerkin method and analyzed the stability of the dual problem.
In Table 2, some numerical methods that have been used to study onedimensional Rosenau equations are summarized. These methods include finite difference and finite element techniques for the space discretization, linear and iterative implicit time discretizations, and time integration with explicit RungeKutta procedures of third and fifthorder accuracy.
 
Equation acronym (E.A.): R = Rosenau, RLW = regularized longwave, B = Burgers, and KdV = Kortewegde Vries. Numerical method: dGM = discontinuous Galerkin method, FDM = finite difference method, FEM = finite element method, QBSPCM = quintic Bsplines collocation method, SOR = successive overrelaxation, CN = CrankNicolson, RK3 = thirdorder accurate RungeKutta method, IL2TL = implicit linear twotime level, IL3TL = implicit linear threetime level, and EB = Euler's backward time discretization. Unless stated otherwise θ = 1. 
Most of the methods presented in Table 2 are secondorder accurate in both space and time, while the one presented in this paper is fourthorder accurate in space and secondorder accurate in time. Moreover, the accuracy of most of the methods presented in Table 2 has been assessed by what is referred to as the method of manufactured solutions; that is, the coefficients of (3) and (4) are found so that the solution of (3) is of a specified travellingwave, for example, hyperbolic secant, type and the initial condition corresponds to the exact solution at . For such an initial condition, is very smooth and does not exhibit steep gradients; as a consequence, a relative small number of grid points are required to obtain very accurate numerical solutions. However, if the initial profile does not correspond to an exact solution of (3) and is negative, the leading part of the initial condition steepens due to the nonlinear advective term and, in the absence of dispersion and diffusion, would result in the formation of a shock wave [16, 17]. For the same type of initial conditions but with dissipation and no dispersion, the initial steepening is eventually balanced by viscous dissipation and a travelling wave which may be referred to as Taylor’s wave is formed [17–19]; such a wave has a finite thickness. On the other hand, for the same initial conditions as the ones discussed above but with dispersion and no diffusion, the initial steepening of the leading part of the initial profile caused by the nonlinear advective term is eventually somewhat balanced by dispersion and a dispersive shock wave or conservative undular bore may form [17, 20, 21].
In this paper, a generalized viscous Rosenau equation that includes linear and nonlinear advective terms is studied numerically by means of a secondorder accurate, linearized CrankNicolson method and threepoint, fourthorder accurate, compact operator discretizations for the first, second, and fourthorder spatial derivatives. The breakup of the initial condition is studied as a function of the linear and nonlinear convective terms and the viscous and two dispersive terms that appear in the equation.
The paper has been arranged as follows. In the next section, the generalized Rosenau equation considered in the study reported here is presented and the linear stability of its linear counterpart is analyzed. This is followed by a section where the numerical method employed to solve the equation and its linear stability are considered. The fourth section presents an exhaustive numerical study of the effects of the linear and nonlinear advective, viscous, and dispersive terms and initial conditions on wave breakup, generation, and propagation. Finally, a short concluding section summarizes the most important findings reported in the paper.
2. Governing Equation
In this paper, the following generalized Rosenau equation is considered:where , , , , and are constants that are associated with the linear and nonlinear advective, dissipative, and third and fifthorder dispersive terms, respectively. Equation (6) is a particular case of (3) and (4), includes both linear and nonlinear terms, and will be solved subject to the following initial and boundary conditions
Before proceeding with the discretization of (6), it is convenient to analyze its linear counterpart.
2.1. Linear Analysis
The linear counterpart of (6) may be written as which, upon introducing , where , and denote the (angular) frequency and wavenumber, respectively, and is the amplitude, becomes which, for real and complex , that is, , where and denote the real and imaginary parts, respectively, of , becomes, upon separating the real and imaginary parts,
The above relations indicate that instabilities arise whenever and this condition requires that for which is not fulfilled if and . Furthermore, (11) and (12) indicate that neither nor is defined if provided that neither nor is zero, respectively.
For small wavenumber, that is, , (11) and (12) indicate that and , thus showing that long waves are stable (the viscosity coefficient ). On the other hand, for very large wavenumbers, that is, very small wavelengths, these equations show that and which indicate that short waves are unstable if . From this analysis, it may be concluded that no linear instabilities occur if and are positive and negative, respectively, and, for these values, is a positive monotonously increasing function of the wavenumber and achieves its minimum value of one for . Moreover, these conditions on and also imply that the phase velocity, that is, , is a monotonously decreasing function of and is equal to , the linear wave speed in (6), for ; on the other hand, for ; that is, there is no damping for .
For , the minimum value of occurs at and is equal to , while the group velocity, that is, , is not defined whenever and and for and , respectively. The group velocity is positive for , nil for , and negative for , where and .
2.2. Invariants
Integration of (6) from to yields where is a constant equal to .
Multiplication of (6) by and integration of the resulting equation from to yield where the total energy density is and contains three semipositive terms if and and, for these values, decreases with time if is not nil. Note that the third term on the righthand side of (15) is positive if .
3. Finite Difference Discretization
Equation (6) may be written as the following system of equations:
Using a secondorder accurate time discretization and linearizing the nonlinear terms with respect to the previous time level, (16) becomes where is the time step, and the superscript denotes the th time level.
Equation (20) together with (17)–(19) represents a linear fourthorder ordinary differential equation for at each time level and is subject to the boundary conditions specified in (8); note that, for example, . However, (20) may not be solved analytically due to the dependence of , , , , , and on , but it becomes a discrete equation upon approximating these variables and , , and by finite differences.
The infinite domain was truncated to a finite one, where and were chosen so that the locations of these boundaries do not affect the numerical results, and the finite domain was discretized into equally spaced nonoverlapping intervals, so that , , , , and is the grid spacing.
In this paper, the following threepoint, compact, fourthorder accurate, finite difference approximations for (17)–(19) have been used:
From (8), and . Furthermore, by assuming that (6) is valid at the boundaries and making use of (8), it is easy to obtain that where the boundary conditions have been applied at the boundaries of the truncated domain. This approximation is valid provided that and are sufficiently far away.
Equation (25) clearly indicates that if , that is, , then , that is, ; these conditions have been used in (22) and (24) to determine and , with , imply that and, together with (8), result in eight boundary conditions for a fourthorder partial differential equation.
It has been found that provided that the solution of (6) is sufficiently smooth near the boundaries of the truncated domain, this overspecification of the boundary conditions does not affect the numerical solution of that equation. A similar behavior was also observed in numerical experiments where and were determined by means of fourthorder accurate, onesided finite difference approximations that make use of five grid points, one located at the boundary and four interior ones, expanding the values of, say, with , in Taylor’s series expansions about and eliminating , , and at from these expansions; an analogous procedure was also used to determine a fourthorder accurate approximation to at the boundaries, except that, in this case, the Taylor series expansions of with were performed with respect to , say, at and , , , , and at were eliminated from these expansions. These onesided approximations, however, result in nonsymmetric, nontridiagonal matrices for and (cf. (22) and (24)) and were found to be less robust than those obtained with (22)–(24) and , especially when the solution of (6) approached the boundaries.
The finite difference method presented in this section is linear at each time step and second and fourthorder accurate in time and space, respectively, and its stencil consists of only three grid points. Moreover, (20) at with (22)–(24) may be written as a block tridiagonal matrix for with , where the superscript denotes transpose that may be solved by means of the block tridiagonal matrix method. Alternatively, one may easily determine , , and from (22)–(24), respectively, where the matrices , , and can be easily determined from those equations and, for example, ; the values of , , and thus obtained can then be substituted into (20) applied at all the interior grid points to obtain a tridiagonal matrix for .
Although the threepoint compact operator method presented here is formally fourthorder accurate in space, it must be noted that this scheme as well as many higherorder ones may result in oscillations in regions where steep gradients exist or the solution is not properly resolved. These oscillations are a consequence of the fact that compact and higherorder methods result in finite difference equations which have more solutions than those of the continuous problem. The spurious solutions of these methods may not cause numerical problems if their magnitude is smaller than those associated with the main roots of the characteristic polynomial of the finite difference equation.
3.1. Linear Stability of the Finite Difference Discretization
Consider the linear counterpart of (6), that is, (9), which upon time discretization may be written as (cf. (20)) and assume that
Substitution of the above expressions in (22)–(24) yields where .
Substitution of (27) and (29) into (26) yields where and the condition of linear stability, that is, , is satisfied provided that which is indeed the case because and . The method is, however, dispersive and the phase of is given by which for , , and is .
4. Results
In this section, some sample results illustrating the numerical solution to (6) are presented for several values of the parameters that appear in that equation. Unless otherwise stated, the initial condition used in the study is which corresponds to a Gaussian function of amplitude and standard deviation .
The numerical experiments reported here were performed with and , and the numerical method was stopped whenever the influence of the boundaries on the solution was noticeable. This explains why, in some presented graphs, the time axis is smaller than in other ones.
Most of the calculations were performed with , that is, , and and correspond to nearly gridindependent results. When either steep gradients or radiation tails were observed, calculations were performed with smaller time steps and grid spacings in order to ensure almost gridindependent results.
The accuracy of the results was assessed in terms of the discrete and norms of the differences between two solutions corresponding to different time steps and/or grid sizes; they were also assessed in terms of the invariants reported before, but the preservation of invariants was found not to be a good indicator of the accuracy of the method because the invariants represent global quantities that are obtained from a spatial integration of the solution at each time level, and numerical quadrature is a smoothing operation.
In many of the results presented here, it was observed that was less than for and , when the invariants were evaluated by means of a secondorder accurate trapezoidal rule (cf. (13)). Further comments on the first invariant are made at the end of this section when assessing the effects of the initial conditions on wave formation.
Before proceeding with the presentation and discussion of results, it is convenient to analyze the different terms that appear in (6) even though the superposition principle is not applicable to such a nonlinear equation. By considering only the first and second terms of (6), that is, , one obtains the wave solution that indicates that remains constant along the characteristic lines , where is a constant and denotes a function. If only the first and third terms of (6) are considered, that is, , the solution of the resulting equation may be written as provided that shock waves are not formed; since the speed of the waves is equal to and, therefore, increases with , wave steepening and shock formation may occur if .
If only the first and fourth terms of (6) are considered, that is, , the solution of the resulting equation may be written as , where and , and decays in both space and time. On the other hand, if only the first and fifth terms of (6) are accounted for, that is, , an integration of this equation yields , the integration of which yields (for ) where and are integration constants; if , the solution to contains trigonometric terms and is not reported here. Note that the homogeneous part of (33) increases as increases unless .
Finally, if only the first and sixth terms of (6) are considered, that is, , an integration of this equation yields , the integration of which (for ) may be written as where , , , and are integration constants and , and indicates that increases as increases unless . The particular solution may be found without much difficulty by means of Lagrange’s variation of parameters, but it is not reported here.
For , the homogeneous solution of is where .
A comparison between (33) and (34) indicates that the former grows faster than the latter as increases if (with and ), that is, if , whereas the opposite holds if . On the other hand, a comparison between (33) and (35) shows that the former grows faster than the latter if (with and ), that is, if , whereas the opposite holds if .
Since, according to the analysis performed in Section 2, for (linear) stability, the contribution of the thirdorder derivative term of (9) is much more important than that of the fifthorder derivative one if ; moreover, since, the shortest wave that can be represented in a grid is and the largest wavenumber that may be resolved is, therefore, , the condition requires that . This implies that must be much larger than . For , a similar condition holds.
Effect of . Figures 1 and 2 show the solution of (6) for two different values of the viscosity coefficient . Figure 1 indicates that the initial Gaussian profile undergoes a large change initially, and its leading edge steepens whereas its amplitude decreases as time increases. The wave profile also widens as time increases; this widening is a consequence of the value of considered in Figure 1, while the steepening of the leading front is a consequence of the nonlinear advective terms, as discussed at the beginning of this section. Figure 1 also shows that, for , .
It should be pointed out that the combs/peaks observed in the amplitude of the leading wave in some figures presented in this paper are due to the fact that only a limited number of grid points and time levels are represented.
As the value of is decreased, the effects of viscosity become less important than the advective and dispersive ones, the wave’s leading front steepening increases, and sawtooth waves may form behind the leading front; these waves have a typical structure that has been observed, for example, in nonlinear acoustics [18], and correspond to dispersive shock waves [17, 20, 21]. The number and amplitude of these sawtooth waves increase as is decreased, and a radiation tail may form as indicated in Figure 2 which corresponds to and the values of the parameters shown in Table 3. Figure 2 also shows the presence of four waves whose amplitude decreases slowly with time.

A similar structure to that of Figure 2 has also been observed for , thus indicating that the number of waves increases from 1 to 4 as is decreased from unity. The width of these waves increases whereas their height decreases as the viscosity coefficient is increased. Figure 2 also shows that the waves formed from the breakup of the initial Gaussian condition used in this study exhibit curved trajectories. Solitary waves of the EW, RLW, and GRLW have also been found to follow curved trajectories in viscous media [45–47].
In Figure 3, the profile is shown at and in order to illustrate the initial adjustment of the wave profile from its initial Gaussian shape, the formation of  or dispersive shock waves behind the leading front, and the number of waves that are formed as functions of the viscosity coefficient. In particular, Figure 3 shows that the amplitude of the leading wave that emerges first from the initial Gaussian profile has a larger amplitude than those that emerge later.
(a)
(b)
Effect of . Figures 4 and 5 illustrate the effects of the linear advection term on the solution of (6). These figures show that there is an initial adjustment of the profile that generates a steep leading front that becomes a solitary wave at later times. The profile behind the trailing edge of the leading wave also steepens and results in the formation of another wave that travels at a smaller speed than the leading one. Behind the second leading wave, a third wave is formed and this process continues as indicated in Figures 4 and 5 until the nonlinear advective terms are smaller than the linear ones. The amplitude of the leading wave is larger than those of the trailing ones; the amplitude of the latter decreases as the distance to the upstream boundary decreases.
Figure 4 also shows (cf. Figures 1 and 2) that the trailing waves are initially of the  or dispersive shock wave type and are caused by the steepening of the profile associated with the nonlinear advective term and dispersion. This is more clearly visible in Figure 5 that corresponds to a linear convective field with .
As is increased, the speed of the leading wave front increases but the number of waves that emerge from the breakup of the initial Gaussian profile decreases; this effect is due to the pushing effect associated with the linear advective term. Although not shown here, similar results to those presented in Figures 4 and 5 have been observed for and the same values of the parameters as those of Figures 4 and 5.
Although not shown here, numerical experiments performed with the same values of the parameters as those of Table 3 but with show similar behavior to that observed in Figures 4 and 5, except that the leading wave speed increases and the leading wave’s front becomes steeper as is increased. This is illustrated in Figure 6 that shows the profile at and 20 for positive and negative values of and indicates that the amplitude of the leading wave is greater than those of the trailing ones. Figure 6 also shows that as the magnitude of for is increased, it takes a longer time to break up the initial Gaussian profile due to the opposing linear convective term.
(a)
(b)