Abstract
The “good” Boussinesq equation is transformed into a first order differential system. A fourth order finite difference scheme is derived for this system. The resulting scheme is analyzed for accuracy and stability. Newton’s method and linearization techniques are used to solve the resulting nonlinear system. The exact solution and the conserved quantity are used to assess the accuracy and the efficiency of the derived method. Head-on and overtaking interactions of two solitons are also considered. The numerical results reveal the good performance of the derived method.
1. Introduction
In recent years, remarkable developments have taken place in the study of nonlinear evolutionary partial differential equations. It is realized that many such equations possess special solutions in the form of pulses which retain their shapes and velocities after interacting with each other. Such solutions are called solitons. Many equations admitting soliton solutions are as follows: sine Gordon and double sine Gordon equations, Schrodinger equation, and KdV, MKdV, and complex modified KdV equations; many research works have been done on these equations [1–8]. Most of the current research is directed to solve coupled nonlinear systems analytically and numerically [9–21]. Solitons are of great interest in many physical areas, as, for example, in dislocation theory of crystals, plasma and fluid dynamics, magnetohydrodynamics, laser and fiber optics, and the study of the water waves.
In this work we will study numerically the “good” Boussinesq (GB) nonlinear equation which provides a balance between dispersion and nonlinearity, that leads to the existence of soliton solutions, similar to the Korteweg-de Vries (KdV) and cubic nonlinear Schrödinger equation [1, 8, 22].
The initial displacement associated with the partial differential equation given in (1) is assumed to take the form [23–25] and the boundary conditions
The solitary wave solution of (1) is given by [3, 14] where is amplitude of the pulse and is an arbitrary parameter. In addition, the double soliton solution [3, 14, 24, 26] can be written as where is the velocity of the th soliton, is the corresponding amplitude, is an arbitrary parameter, and is the initial position. The exact solution (5) represents two solitary waves located initially at the positions and moving to the right or left according to their velocities.
Many research works on Boussinesq equation have been developed. Analytical solution of this equation was studied by many authors, such as the construction of -soliton solutions using the bilinear form [4], multiple soliton solutions for the GB equation using a simplified version of Hirota method [27] and Adomian decomposition method [28]. Construction of soliton solutions and periodic solution of Boussinesq equation by modified decomposition method are given in [29, 30]. A variational iteration method is developed for GB equation [31]. A solitary wave solution of the Boussinesq equation with power law nonlinearity is derived in [32].
Many numerical methods have been developed for solving the Boussinesq equation, such as Petrov-Galerkin method [19]. Bratsos et al. [23, 24, 33] have developed finite difference schemes and have considered the Bad Boussinesq (BB) and GB equations. El-Zhoheiry has developed an implicit finite difference scheme [34]. Aydin and Karasözen [35] constructed second order symplectic and multisymplectic integrators for the GB equation using the two-stage Lobatto IIIA-IIIB partitioned Runge-Kutta method. Daripa and Hua [36] developed finite difference method for BB equation of second order in time and space, where they have used the filtration and regularization techniques to control the growth of the errors arising from the instabilities and provide better approximate solutions of this equation. Ismail and Bratsos [25] have derived a predictor corrector scheme for the GB and BB equations; the scheme is fourth order in time and second order in space, and it is conditionally stable. Matsuo [37] also has derived conservative finite difference schemes for certain classes of nonlinear wave equations with some applications for the nonlinear Klein-Gordon and GB equations. Mohebbi and Asgari [26] also have solved the GB equation using a fourth order time stepping schemes with combination of discrete Fourier transform. Split step Fourier method is also used to solve Boussinesq-type equations. Dehghan and Salehi [38] have derived a mesh free method for the classical Boussinesq equation.
The paper is laid out as follows. In Section 2, the GB equation is transformed into a first order differential system in time; a finite difference scheme is derived for this system. The method is analyzed for accuracy and stability. In Section 3, a linearization technique is used to convert the nonlinear system obtained into a linear block tridiagonal system. Numerical tests are presented in Section 4. Concluding remarks are given in Section 5.
2. The Numerical Method
In order to derive a highly accurate method, we transform the GB equation (1) into the first order differential system in time as [19] with initial conditions and boundary conditions
By using the boundary conditions (10), the system (7) and (8) has the conserved quantity [24–26, 38]
In order to derive a numerical method for solving the system (7) and (8), the region with its boundary consisting of the lines , and the axis is covered with a rectangular mesh of points with coordinates [25] where and are the space and time increments, respectively. We denote the theoretical solution of (7) and (8) at the typical mesh point by and , and the numerical solution by and . The fourth finite difference formula [39] is used to approximate the second order space derivative which appears in our system. By making use of (13) and the implicit midpoint rule, the finite difference nonlinear implicit scheme is obtained. The schemes (15) and (16) can be written as where . On expansion of the central differences using (14), this will lead us to the difference scheme By imposing the boundary conditions, a coupled nonlinear system is obtained, which can be solved by any iterative scheme. This system can be solved by Newton's method or fixed point method. Linearization method can also be used.
2.1. Accuracy of the Scheme
To study the accuracy of the proposed schemes (19) and (20), the numerical solutions and are replaced by the exact solutions and , respectively. By making use of the following expansions about the point : into (19) to obtain
By making use of the differential system (7) and (8), all terms inside the brackets in (22) are zero. Similar analysis can be done for (20). This will lead us to the conclusion that the derived scheme is of second order in time and fourth order in space, that is, . The numerical results confirm this.
2.2. Stability of the Scheme
To study the stability of the derived scheme, von Neumann stability analysis will be used. We consider the linear version of the proposed schemes (17) and (18) which can be given as where is a constant quantity defined by
To apply von Neumann method, we assume that the solutions of (22) and (23) are where , , and are constants. By substituting these solutions into (23) and (24), and after some manipulation, we get the system where . Equations (27) can be written in a matrix vector form as where
The von Neumann necessary condition for the stability of this system is the maximum eigenvalue of value of the amplification matrix in (28) is to be less than or equal to one. By direct calculation, the eigenvalues of the matrix are with modulus equal to one, and hence the proposed scheme is unconditionally stable.
3. Linearization Technique
In order to avoid the solution of the nonlinear systems (19) and (20), a linearization technique will be developed to overcome this difficulty. Using Taylor's series expansion of the nonlinear term about , we obtain and this will lead us to the following approximation: which preserves second order temporal accuracy. By replacing the exact solution by the approximating one, we get
By substituting (33) into (18), we obtain the linearly implicit finite difference scheme
By expanding the central difference operator in (34), we obtain the linearized scheme
By utilizing the boundary conditions, the finite difference scheme (35) produces a linear block tridiagonal system in the unknowns , which can be solved directly by Crout's method. The obtained difference scheme is still of second order in time and fourth order method in space, and it is unconditionally stable.
4. Numerical Results
In all experiments we use the following values: , , , , and . We study the accuracy of the proposed method by calculating the infinity error norm The numerical solution and trapezoidal rule are used to calculate the error and the conserved quantity [19, 25, 26, 37]. All numerical results in this section are obtained from the solution of the nonlinear schemes (19) and (20) using Newton's method. Linearization method is used for comparison purpose only.
4.1. Single Soliton
To test the derived method, we consider the initial condition where (37) represent soliton and kink solution, respectively. The parameters , , and are used in this test. The numerical results for the nonlinear and linearized schemes are given in Tables 1 and 2, respectively. The numerical results (amplitude and the conserved quantity) obtained display the high accuracy of the proposed method. The execution time required for the nonlinear scheme is sec. compared to sec. for the linearized scheme. The numerical results produced in Table 1 are more accurate than the one in Table 2, and this is due to the approximation of the nonlinear term in the linearization process. A comparison of some existing methods is given in Table 3, which indicates that our method is the most accurate one. In Figures 1 and 2, we display the numerical solutions and moving to the right for .
To test whether the proposed numerical scheme exhibits the expected convergence rate in space, we perform some numerical experiments for various values of and fixed value of . In these experiments we choose to ensure that the temporal error is negligible. The rate of convergence for the scheme is calculated using the formula where is the infinity error norm. The errors and the rate of convergence for are given in Table 4.
4.2. Head-On Soliton Interaction
To study the head-on collision of two solitons, we choose the initial conditions where
In [19] it was reported that solution blew up numerically when . By choosing the parameters , , and . The initial conditions represent two solitons and two kinks located at and . In Figures 3 and 4, we display the interaction scenario for to ; we can easily observe that the two solitons and the two kinks have been separated completely after the interaction and restored their original shapes and velocities. The velocity and the conserved quantity are given in Table 5. It is noted that when the two waves overlap (Figures 3 and 4), the joint amplitude is greater than the sum of individual amplitudes; this is in full agreement with [37]. By choosing , the situation dramatically changes; both Figures 5 and 6 end around ~. Newton's method fails to find the solution there. This agrees with the blowup results in [19, 37].
4.3. Overtaking Soliton Interaction
In this test we will study the interaction of two solitons moving in the same direction. In this case we choose the initial conditions where
The following parameters are selected:
The initial conditions represent two solitons and two kinks are initially located at and , moving in the same direction with velocities and . The usual interaction has taken place and the faster wave interacted and separated from the slower one and left it behind; the interaction scenario is given in Figures 7 and 8. In Table 6, we track the amplitude and the conserved quantity during the interaction scenario. In Figure 9, we display the contours of the numerical solution . Many numerical tests have been conducted for different amplitudes. We have noticed that the interaction occurred for , , which is completely different from the head-on interaction case, where the solution blew up numerically when , . This type of interaction for the GB equation seems to be reported in the literature for the first time as far as the authors know. We have noticed that when the two waves overlap (Figures 7 and 8), the joint amplitude is less than the sum of individual amplitudes.
4.4. Birth of Solitons
To study the evolution of arbitrary pulse [19], we choose the initial conditions
This test is carried out for values of , mainly in the vicinity of , the theoretical value of the amplitude for . We considered two cases and . In the first case we choose ; we have noticed that as time evolves, the initial pulse split into two solitons moving in the opposite directions with equal amplitudes and velocities ; see Table 7. In Figures 10 and 11, we display the splitting scenario and some dispersive oscillations which occurred after splitting. In the second case we choose ; as time evolves, the initial pulse is split into two solitons moving in the opposite direction with amplitudes ; see Table 8. Figures 12 and 13 display the splitting scenario; no dispersive oscillation is observed in this case (linearized scheme failed to produce the numerical solution in this case). Many numerical tests have been conducted for different amplitudes; we have noticed that the numerical solution blew up if .
5. Concluding Remarks
In this paper, a highly accurate finite difference method to solve GB equation is proposed. The main idea is to transform the GB equation into a first order differential system in time. The method is fourth order in space and second order in time, and it is unconditionally stable. The method produced a coupled nonlinear system; Newton's and linearization methods are used to solve it. The nonlinear scheme produced more accurate results than the linearized scheme. Numerical results for single soliton, head-on, and overtaking solitons interactions are given; the most interesting phenomena we observed were when two solitons collide; if both solitons are small enough, they pass through each other like the other usual solitons do, but when they exceed some limit, the solution blows up at the collision, even if both amplitudes are smaller than for being stable solitons [14, 15]. The overtaking solitons interaction is also discussed as a new numerical test. Elastic interaction occurred for large amplitudes (, , and this is contrary to the head-on interaction, where blowup is observed for moderate values of the amplitudes (, . Birth of solitons is also observed for initial pulse with amplitude . A blowup was observed when and this is in full agreement with the existing methods [14, 15].
Conflict of Interests
The authors declare that they have no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors would like to express their sincere thanks to the referees for their valuable comments on the improvement of the paper.