Abstract
A Galerkin method for a modified regularized long wave equation is studied using finite elements in space, the Crank-Nicolson scheme, and the Runge-Kutta scheme in time. In addition, an extrapolation technique is used to transform a nonlinear system into a linear system in order to improve the time accuracy of this method. A Fourier stability analysis for the method is shown to be marginally stable. Three invariants of motion are investigated. Numerical experiments are presented to check the theoretical study of this method.
1. Introduction
Much physical phenomena are described by nonlinear partial differential equations. Most of these equations do not have an analytical solution, or it is extremely difficult and expensive to compute their analytical solutions. Hence a numerical study of these nonlinear partial differential equations is important in practice. The regularized long wave (RLW) equation where is a positive constant, is a nonlinear evolution equation, which was originally introduced by Peregrine [1] in describing the behavior of an undular bore and studied later by Benjamin et al. [2]. This equation plays an important role in describing physical phenomena in various disciplines, such as the nonlinear transverse waves in shallow water, ion-acoustic waves in plasma, magnetohydrodynamics waves in plasma, longitudinal dispersive waves in elastic rods, and pressure waves in liquid’s gas bubbles. Many numerical methods for the RLW equation have been proposed, such as finite difference methods [3, 4], the Galerkin finite element method [5–8], the least squares method [9–11], various collocation methods with quadratic B-splines [12], cubic B-splines [13] and septic splines [14], meshfree method [15, 16], and an explicit multistep method [17].
The RLW equation is a special case of the generalized regularized long wave (GRLW) equation where and are positive constants and is a positive integer. Some numerical methods [18–24] for the GRLW equation have been also presented, such as a finite difference method [18], a decomposition method [20], and a sinc-collocation method [23]. Another special case of the GRLW equation is called the modified regularized long wave (MRLW) equation in which [25]. Some authors have studied the MRLW equation using various numerical methods, such as a finite difference method [26], a collocation method [27], a spline method [28, 29], and the Adomian decomposition method [30].
In this paper, we study a Galerkin method for the MRLW equation by using linear finite elements in space and extrapolation to remove the nonlinear term. We discuss the properties of this method and compare its accuracy with previous studies. The interaction of two and three solitons is also studied. Moreover, the propagation of the Maxwellian initial condition is simulated.
The outline of this paper is as follows. In the next section, the governing equation, its analytical solution, and three invariants are given. In Section 3, we propose the numerical method, including semidiscretization in space and full discretization in space and time. In Section 4, we present a stability analysis. In Section 5, some numerical experiments are presented. Finally, we give a brief conclusion in Section 6.
2. Governing Equation
We assume to be an open, bounded subset of and to be a final time. We set , and consider the following MRLW equation: where is the initial datum, is a positive constant, and the subscripts and denote differentiation in space and time, respectively. The physical boundary condition requires as .
The MRLW equation (3) has the exact solution [25]: where and and are arbitrary constants. Furthermore, (3) possesses three invariants of motion corresponding to conservation of mass, momentum, and energy [25]: These invariants are used to check the conservative properties of a numerical method, especially for problems without an analytical solution and during collision of solitons.
3. Numerical Methods
Set . Applying Green’s formula to problem (3) and using the boundary condition in the definition of , we derive the variational form of problem (3). Find such that where
3.1. Semidiscretization in Space
Let us consider a uniform mesh with the mesh size , , which consists of points . Then we obtain We make the transformation , , , in order to transfer an element into a standard interval. From (8), we derive the following equation:
Now, we define the finite dimensional subspace , , where are the linear basis functions on each element. Then the semidiscrete scheme for problem (3) is formulated as follows. Find such that
The variation of over the element , is expressed as
For , we choose in problem (11) and substitute (12) into (11). Then an element’s contribution is obtained in the form of where the symbol () denotes differentiation with respect to , which, in matrix form, is given by where are relevant nodal parameters. The element matrices are
For the element , the indices and take only the values and so that the matrices , , , and are 2 × 2:
Assembling contributions from all elements leads to the following matrix form of the coupled nonlinear ordinary differential equations: where contains all the nodal parameters. The four assembled matrices are tridiagonal. The general row for each matrix has the following form:
3.2. Full Discretization
We now consider a fully discrete scheme. Let be a uniform partition of into subintervals , , with length . For a generic function of time, set . We use the Crank-Nicolson discretization method in (17) with Then (17) can be written as with . This equation is symmetric about the point , and one should, therefore, expect second order accuracy in time. However, the linearization decreases the order of the time discretization error to . This drawback can be overcome by using an extrapolation technique in the linearization of the nonlinear coefficient . We choose , for , and define Now, the Crank-Nicolson method can be shown to produce an error of order . Combining with the error of finite element method, the Crank-Nicolson Galerkin finite element method of (6) produces an error of order overall [17].
We can obtain the following recurrence relation at point : where and are given by
As (22) is used only for , we must compute by another method. For this, we choose a predictor-corrector method. Let in (21); we can get the first approximation with replaced by . Then we use to substitute in (21) and regard the result as the final approximation . Or, more precisely, the concrete procedure is as follows. First, set then we calculate finally, we have
Thus we can use the Thomas algorithm to solve the linear algebraic equations.
In order to improve the accuracy in time, we also use Runge-Kutta discretization method. The fourth order Runge-Kutta discretization method in (17) can be described as follows:
4. Stability Analysis
A linear stability analysis is made of the growth factor of the error in a typical Fourier mode of amplitude : where is a mode number and is the element size. The von Neumann stability method can only be applied to a linear scheme, so we linearize the nonlinear term by assuming that in this term is locally constant. Substituting (28) into (22) gives where Taking the modulus of (30) gives Therefore, our method is marginally stable.
5. Numerical Experiments
In this section, we present some numerical tests to check the efficiency and accuracy of our method, which are the propagation of single soliton and collision of two and three solitons at different time levels. Finally, we investigate the development of the Maxwellian initial condition into solitary waves.
To illustrate the accuracy of the present method, we use - and -error norms to compare the numerical solution with the exact solution, which show the mean and maximum differences between the numerical and analytical solutions. The quantities , , and measure the conservation laws of our method during propagation.
5.1. Single Soliton
The analytical values of the three variants are
For the purpose of comparison with the previous work [28], we choose the parameters , , , , and . Then the amplitude is and the analytical values of three invariants are , , and . The relevant numerical results which applied the Crank-Nicolson scheme in (22) and the Runge-Kutta scheme in (27) are presented in Tables 1 and 2, respectively. The profiles of solitary waves at time = 5 and time = 10 are given in Figure 1 and the error distribution in the wave profile in the three-dimensional view at different time level is illustrated in Figure 2. Table 3 presents the variants, norm, and norm at time = 10 against the quadratic B-splines and the cubic B-splines in [28].
We also choose the parameters , , , and , with the range . Then the amplitude is and the analytical values of three invariants are , , and . Tables 4 and 5 show the results about the three invariants, the norm, and the norm at different times and compare these results with those by the collocation method with cubic B-splines in [27] using the Crank-Nicolson scheme in (22) and the Runge-Kutta scheme in (27), respectively. The changes of the variants are satisfactorily small, though our result for this case is not more accurate than the result in [27], where the simulation is done up to . The profile of the approximate solution is given in Figures 3 and 4 shows the error profile of the single solitary waves.
The Runge-Kutta Galerkin method algorithm has been run for different space steps with a fixed time step . The error and convergence order of space for the - and -norm are recorded in Tables 6 and 7, respectively, which verify that the - and -error can achieve the order in space. And the algorithm has also been run for different time steps with a fixed space step . The error and convergence order of time for the - and -norm are recorded in Tables 8 and 9, respectively.
5.2. Collision of Two Solitons
In this section, we study the interaction of two solitary waves with the initial conditions given by a linear sum of two separate solitary waves of various amplitudes: where , , , and and are arbitrary constants. The analytical values for the conservation laws in this case have the following form: In our numerical scheme, we choose , , , , , and , with interval . Then the amplitudes are in ratio 2 : 1, where and the analytical values for the conservation are , , and . The changes of the invariants are satisfactorily small, since the changes of the invariants , , and are 1.3295 × 10−4, 2.2259 × 10−5, and 6.8989 × 10−5, respectively. The results are recorded in Table 10. A graph of the two solitary waves collision, plotted in a three-dimensional view at some discrete times, is shown in Figure 5.
5.3. Collision of Three Solitons
In this section, we study the interaction of three solitary waves with the initial conditions given by a linear sum of three separate solitary waves of various amplitudes: where , , , and and are arbitrary constants. The analytical values for the conservation laws in this case have the following form: In this case, we choose paraments as , , , , , , , and , with interval . Then the amplitudes are in ratio 4 : 2 : 1, where and the analytical values for the conservation are , , and . The changes of the invariants from the initial variants approach zero throughout and agree with the analytical values for the three invariants as presented in Table 11, which indicates that our scheme is satisfactorily conservative. The geometry of the initial state and the profiles at time , , and is shown graphically in Figure 6, respectively.
(a) Time = 0
(b) Time = 10
(c) Time = 20
(d) Time = 35
5.4. The Maxwellian Initial Condition
Finally, we consider the development of the Maxwellian initial condition into a train of solitary waves. We choose different values , , and in our numerical scheme. The comparison of the three variants , , and with the earlier result in [29] is presented in Table 12, here Crank-Nicolson method is used for time discretization. For , the changes of the variants , , and with respect to the initial values are 3.9492 × 10−6, 6.9575 × 10−4, and 3.1509 × 10−3, respectively, whereas they are 2.8210 × 10−6, 2.60000 × 10−4, and 7.0771 × 10−4 in [29]. For , the variants are changed by 7.8982 × 10−6, 6.8214 × 10−5, and 6.1410 × 10−4 in this case. The development of the Maxwellian initial condition is shown in Figure 7 with different parameter values, respectively. The smaller there is, the more the number of solitary waves will form. The simulations are done up to time = 10 in this case. From the previous work [25], the total number of solitary waves and the values have the approximate relation
| (a) |
| (b) |
| (c) |
6. Conclusion
In this paper, we have developed a Galerkin linear finite element method to investigate the propagation of solitons and their interactions governed by the nonlinear MRLW equation. An extrapolation technique has been used to improve the accuracy in time for this numerical method. A linear stability analysis shows that this method is marginally stable. The high efficiency and accuracy of our method is tested by the numerical experiments: propagation of single soliton, collision of two and three solitons, and development of the Maxwellian initial condition into solitary waves. Moreover, this method numerically satisfies the conservation laws of mass, momentum, and energy.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The work was supported by the National Natural Science Fund of China (11371289) and the China Scholarship Council (201206285005).