In this work, our main purpose is to develop of a sufficiently robust, accurate and efficient numerical scheme for the solution of the regularized long wave (RLW) equation, an important partial differential equation with quadratic nonlinearity, describing a large number of physical phenomena. The crucial idea is based on the discretization of the RLW equation with the aid of a combination of the discontinuous Galerkin method for the space semi-discretization and the backward difference formula for the time discretization. Furthermore, a suitable linearization preserves a linear algebraic problem at each time level. We present error analysis of the proposed scheme for the case of nonsymmetric discretization of the dispersive term. The appended numerical experiments confirm theoretical results and investigate the conservative properties of the RLW equation related to mass, momentum and energy. Both procedures illustrate the potency of the scheme consequently.
PACS Codes: 02.70.Dh, 02.60 Cb, 02.60.Lj, 03.65.Pm, 02.30.-f.
MSC: 65M60, 65M15, 65M12, 65L06, 35Q53.
Keywords:discontinuous Galerkin method; regularized long wave equation; backward Euler method; linearization; semi-implicit scheme; a priori error estimates; solitary and periodic wave solutions; experimental order of convergence
We are concerned with a proposal of a sufficiently robust, accurate and efficient numerical method for the solution of scalar nonlinear partial differential equations. As a model problem, we consider a regularized long wave (RLW) equation firstly introduced by Peregrine (in ) to provide an alternative description of nonlinear dispersive waves to the Korteweg-de Vries (KdV) equation. As a consequence of this, the RLW can be observed as a special class of a family of KdV equations.
The RLW equation contains a quadratic nonlinearity and exhibits pulse-like solitary wave solutions or periodic waves; see . It governs various physical phenomena in disciplines such as nonlinear transverse waves in shallow water, ion-acoustic waves in plasma or magnetohydrodynamics waves in plasma. Since the RLW equation can be solved by analytical means in special cases, the proposed numerical methods can be easily verified. Several numerical studies of the RLW equation and its modified variant have been introduced in the literature, from finite difference methods , over collocation methods [4,5], to finite element approaches [6,7], or Galerkin methods , and in references cited therein.
In this paper, we present a semi-implicit scheme for the numerical solution of the RLW equation based on an alternative approach to the commonly used methods. The discontinuous Galerkin (DG) methods have become a very popular numerical technique for the solution of nonlinear problems. DG space semi-discretization uses higher-order piecewise polynomial discontinuous approximation on arbitrary meshes; for a survey, see [9,10] and . Among several variants of DG methods, we deal with the nonsymmetric variant interior penalty Galerkin discretizations; see . The discretization in time coordinate is performed with the aid of linearization and the backward Euler method, sidetracking the time step restriction well known from the explicit schemes, proposed in . Consequently, we extend the results from , and the attention is paid to the a priori error analysis of the method with the aid of standard techniques introduced in  and .
The rest of the paper is organized as follows. The problem formulation and its variational reformulation are given in Section 2. Discretization, including space semi-discretization and fully time space discretization, is considered in Section 3. The Section 4 is devoted to a priori error analysis. Finally, in Section 5, the theoretical results are illustrated by numerical tests on a propagation of a single solitary wave and experimental orders of convergence are computed for piecewise linear approximations together with invariant quantities of the RLW equation.
2 Regularized long wave equation
Let be a bounded open interval and be a final time. We set and consider the following initial boundary value problem: Find such that, for all ,
where constant parameters and are related to the amplitude of the wave and long-wavelength, respectively. From the mathematical point of view, problem (1)-(3) represents the regularized long wave equation equipped with a set of two generally nonhomogeneous Dirichlet boundary conditions (2) and with the initial condition . These given data have to satisfy the compatibility conditions prescribed at both endpoints of the domain Ω, i.e.,
The whole system (1)-(4) was found to have single solitary or periodic traveling wave solutions; for details, see .
Remark 1 In the case of a single solitary wave propagation, the homogeneous Dirichlet boundary conditions in (2) arise from the asymptotic behavior of the exact solution u, and the endpoints a and b are chosen large enough so that the boundaries do not affect the single solitary wave during its propagation up to final time T.
In what follows we use the standard notation for function spaces and their norms and seminorms . Let be an integer and . We use the well-known Lebesgue and Sobolev spaces , , Bochner spaces of functions defined in with values in a Banach space X and the spaces of k-times continuously differentiable mappings of the interval [0, T] with values in X. By we denote the subspace of all functions satisfying . To this end, we use the following notation for a scalar product in by
for a norm in by , for a seminorm in by and for a norm in by . It is a known fact that is a norm on equivalent to
A sufficiently regular solution satisfying (1)-(3) pointwise is called a classical solution. Now, we are ready to introduce the concept of weak formulation. Firstly, we recall the definition of a bilinear dispersion form and a nonlinear convection form from , i.e.,
where symbol stands for the function on Ω such that , and function in (7) represents the physical flux.
Definition 1 We say that u is a weak solution of problem (1)-(3) if such that and the following conditions are satisfied:
Remark 2 In order to unify the definition of the weak solution (8), we consider nonhomogeneous Dirichlet boundary conditions instead of the second parallel analysis of periodic-type solutions with the aid of Sobolev spaces of periodic functions with period satisfying .
Further, to carry out the error analysis later, we need to specify additional assumptions on the regularity of a solution of continuous problem (1)-(3). Therefore, we assume that the weak solution u is sufficiently regular, namely
Let ( ) be a family of the partitions of the closure of the domain Ω into N closed mutually disjoint subintervals with length and the symbol stands for an index set . Then we call a partition with a spatial step and interval an element. By we denote the set of all partition nodes of Ω, i.e., . Further, we label by the set of all inner nodes. Obviously, .
We additionally assume that the partitions satisfy the following condition.
Assumption (M) are locally quasi-uniform:
The condition (12) in fact allows to control a level of the mesh refinement if adapted meshes are used.
DG methods can handle different polynomial degrees over elements. Therefore, we assign a local Sobolev index and a local polynomial degree to each . Then we set the vectors and . Over the triangulation , we define the so-called broken Sobolev space corresponding to the vector s as
with the norm
and the seminorm
where and denote the standard norm and the seminorm on the Sobolev space , .
The approximate solution of variational problem (8) is sought in a finite dimensional space of discontinuous piecewise polynomial functions associated with the vector p by
where denotes the space of all polynomials of degree on , .
Let us denote and . Then we can define the jump and average of v at inner points of the domain Ω by
By convention, we also extend the definition of jump and mean value for endpoints of Ω, i.e., , , and . In case that are arguments of or , we usually omit these arguments , and write simply and , respectively.
3.1 DG semi-discrete formulation
Now, we recall the space semi-discrete DG scheme presented in . First, we multiply (1) by a test function , integrate over an element and use integration by parts in the dispersion term and convection term of (1) subsequently. Further, we sum over all and add some artificial terms vanishing for the exact solution such as penalty and stabilization terms, which replace the inter-element discontinuities and guarantee the stability of the resulting numerical scheme, respectively. Consequently, we employ the concept of an upwind numerical flux (see ) for the discretization of the convection term and end up with the following DG formulation for the semi-discrete solution , introduced in  as a system of ordinary differential equations, i.e.,
where forms and stand for the semi-discrete variants of forms (6) and (7), i.e.,
The crucial item of the DG formulation of the model problem is the treatment of the convection part. We proceed analogously as in , where the convection terms are approximated with the aid of the following numerical flux through node in the positive direction (i.e., outer normal is equal to one):
where and the choice of and for boundary points has to satisfy the prescribed Dirichlet boundary conditions; for more details, see .
In what follows, we shall assume that the numerical flux has the following properties.
(H1) is Lipschitz-continuous with respect tou, v:
(H2) is consistent:
(H3) is conservative:
One can see that the numerical flux H given by (21) satisfies conditions (H2) and (H3) and is Lipschitz-continuous on any bounded subset of .
A particular attention should be also paid to the treatment of the dispersion terms, which include an artificially added stabilization in the form , in order to guarantee the stability of the numerical scheme. In our case, where this stabilization is added with a positive sign (+), we speak of the nonsymmetric interior penalty Galerkin method.
In the end, the semi-discrete DG scheme is completed with the weighted penalty
which replaces the inter-element discontinuities and guarantees the fulfillment of the prescribed boundary conditions.
The penalty parameter function for a nonsymmetric variant is defined in spirit of  as
In order to simplify the notation, we introduce the form
which is bilinear due to (19) and (25). Consequently, we can here define the semi-discrete solution of problem (8).
Definition 2 We say that is a semidiscrete solution of problem (8) if and the following conditions are satisfied:
3.2 Semi-implicit linearized DG scheme
In order to obtain the discrete solution, it is necessary to equip the scheme (28) with suitable solvers for the time integration. In , we have proposed a semi-implicit time discretization based on the backward Euler scheme with the linearized convection form which is suitable for avoiding the strong time step restriction of explicit time schemes as well as for the preservation of linear algebraic problems at each time level.
We now partition as , denoting each time step by and letting stand for the approximate solution of , , . The linearization of the physical flux f is treated in spirit of  as
which implies the splitting of a convection form in the following way:
where and represent the corresponding linearized and nonlinear parts of the original numerical flux given by (21); for more details, see . One can easily observe that the form is linear in its second and third argument and the form is in fact an original convection form (20) with half the amount of the physical flux but from the previous time level.
The fully discrete solution of problem (18) via the aforementioned semi-implicit approach is defined in following way.
Definition 3 Let be a partition of the time interval and , . We define the approximate solution of problem (8) as functions , , , satisfying the following conditions:
Discrete problem (33) is equivalent to a system of linear algebraic equations at each time instant . In what follows, we shall be concerned with the analysis of method (33).
Lemma 1Discrete problem (33) has a unique solution.
Proof We rewrite problem (33) in the following way. For , and , we find such that
Using the definitions (27) and (31), one can see that is a bilinear form on the finite dimensional space and is a linear functional. Moreover, the form is coercive, i.e.,
Hence, equation (34) has a unique solution . □
4 A priori error analysis
For error analysis and in experiments, we consider for all . Thus we denote . Now we would like to analyze the error estimates of the approximate solution , , obtained by method (33). For simplicity, we consider a uniform partition , , of the time interval with time step , where is an integer.
Let be the standard -interpolation of , ( ) satisfying (cf.) for all , ,
for a generic constant independent of h and v. We set
Then the error can be expressed as
Setting in (33), we get
On the other hand, from (18) it follows
Multiplying (43) by τ, subtracting from (42) and using again the linearity of the form , we get
Since , we can rewrite equation (44) in the following way:
For the term on the right-hand side of equation (45), we use decomposition and linearity of the form . Together with (30), we finally get
For next estimates, we use the following lemmas.
Lemma 2Under assumptions (R) for , the following hold:
wherecis a generic constant independent ofhandτ.
Proof The proof of these standard estimates can be found, for instance, in . □
Lemma 3Under assumptions (R), (H) and for , the following hold:
wherecis a generic constant independent ofhandτ.
Proof Again, one can find the proof of these estimates in . □
Since is always nonnegative, applying previous lemmas gives us
Multiplying by 2, applying the Young inequality and using the definition of the form , we obtain
If we take into account that , the previous inequality can be rewritten as
where we denoted and
Let us now introduce the so-called energy norm
and the norm
Denoting and from (58), it follows
In order to finish our estimates, we require a fulfillment of the following technical assumption.
(T1) There exists such that .
If assumption (T) is fulfilled, then . Thus we can also reformulate assumption (T) so that .
Thus, let us assume that assumption (T) holds, then
with . Consequently,
and since , we have , i.e.,
where we used a straightforward estimate . We can notice that due to the presence of the factor μ in front of the function on the left-hand side of (65), we lost μ in denominators of .
Now we are ready to formulate the main theorem.
Theorem 1Let assumptions (M), (H), (R) and (T) be satisfied, then there exists a constant such that
where is defined by (60).
Proof Since , the statement of the theorem comes from (65) and the fact that . Then we set
Remark 3 Theorem 1 implies that the error of our method is in both energy and -norm. However, as we will see in the next section, the error estimate in the -norm is suboptimal with respect to h.
Remark 4 The dependency C on μ in the expression (68) (choice of θ depends on μ) can be removed by applying the so-called continuous mathematical induction mentioned in . This is useful namely in the cases when convection terms dominate, i.e., . Consequently, in these cases assumption (T) can be weakened to a CFL-like condition for suitable .
5 Numerical experiments
In this section we shall numerically verify the theoretical a priori error estimates of the proposed semi-implicit method (33) for the cases of propagation of both a single solitary wave and periodic waves.
We verify numerically the convergence of the method in the -norm and the energy norm given by (60) with respect to time step τ and mesh size h. The computational errors are evaluated at certain time instants during all computations in the corresponding norms, i.e.,
where is a prescribed exact solution at time lτ and is the numerical solution at time level lτ obtained by the semi-implicit scheme (33) with constant time step τ on the uniform grid with mesh size h. We suppose that errors behave according to the formula
The constants , , are independent of τ and , , are independent of h. The values , , , are the orders of accuracy of the method in the corresponding considered norms. We define the experimental order of convergence (EOC) by
5.1 Single solitary case
The RLW equation has the following analytical single solitary wave solution given by
which represents a single solitary wave of amplitude 3c, traveling with the velocity in a positive x-direction and located initially at the point . The initial condition is extracted from the exact solution (74) and homogeneous Dirichlet boundary conditions are set.
In order to compare our semi-implicit approach to the schemes given in [6,20] and , we set the parameter values , , . The run of the algorithm is carried out up to time over the problem domain . The resulting linear algebraic problems (33) are solved by the restarted GMRES method. Figure 1 captures the development of approximation solutions of a single solitary wave from an initial condition to the terminal time T for a piecewise linear approximation with time step and mesh size . The similar plots are also obtained for another combination of τ and h as we consider below.
Figure 1. The 3D plot of approximation solutions of a single solitary wave (left) and corresponding isolines in space-time domain (right).
5.1.1 Convergence with respect to h
First, we investigate the convergence of the method with respect to h. In order to restrain the discretization errors with respect to time step τ, we use a sufficiently small time step . Numerical experiments are carried out with the use of piecewise linear ( ) approximations on five consecutive uniformly refined meshes having 125, 250, 500, 1,000 and 2,000 elements.
Tables 1 and 2 show computational errors in the -norm and the energy norm at four time instances , , and , the corresponding EOC during all computations. Since the exact solution is sufficiently regular over Ω, it follows from Remark 3 that the theoretical error estimates are of order . On the other hand, we observe that the numerical experiment of propagation of a single solitary wave indicates a better behavior of EOC in the -norm, which is expected to be asymptotically for piecewise linear ( ) approximations. These observations also correspond with the finite element approach from , where the same example was studied.
Table 1. Single solitary case: Computational errors in the -norm and experimental orders of convergence for approximation on a consequence of meshes at time instancest( )
Table 2. Single solitary case: Computational errors in the energy norm and experimental orders of convergence for approximation on a consequence of meshes at time instancest( )
Further, the results for EOC in the energy norm are in a quite good agreement with derived theoretical estimates; in other words, this technique produces an optimal order of convergence . Finally, both estimates in Theorem 1 confirm the well-know attribute of DG schemes from the class of convection-diffusion problems, cf. and .
5.1.2 Convergence with respect to τ
Secondly, we verify experimentally the convergence of the method in the -norm and the energy norm with respect to time step τ. In order to restrain the discretization errors with respect to h, we use a fine mesh with 2,000 elements with piecewise linear approximation.
The computations were carried out with five different time steps τ, see Table 3. The computational error is evaluated at final time in the -norm and the energy norm, respectively. We observe that both computational errors have EOC of order in the corresponding norms, which is again in a good agreement with derived theoretical results.
Table 3. Single solitary case: Computational errors in the -norm and the energy norm for approximation with respect to time step ( )
5.1.3 Invariant conservation quantities
Similarly as in , we shall monitor the three conservation quantities for the propagation of the single solitary wave corresponding to mass
with respect to the run of the proposed algorithm. The analytical values for the invariants on the entire real domain are given (in ) by
Moreover, for the purpose of a more accurate comparison with reference results, we introduce the discrete -norm defined by
assessing the accuracy of the method by measuring the difference between the numerical and analytic solutions and u, respectively.
Table 4 records the invariant quantities together with errors , computed on the finest space-time grid and compares obtained results with several previously presented schemes given in [6,7,20] and . All three conservation quantities are kept almost constant thus they illustrate the suitability of the proposed scheme for this problem. The obtained satisfactory results correspond to the reference ones from [6,7,20] and . Moreover, we append the recent results from  in order to compare the incomplete variant of stabilization in dispersion terms with the nonsymmetric one.
Table 4. Single solitary case: Computed invariant quantities and errors in the -norm and the -norm ( , )
5.2 Periodic case
The family of periodic solutions of the RLW equation may be analytically written as (cf.)
and the parameters , and given by
where is the Jacobi elliptic function and stands for the elliptic modulus; for definitions and other properties, see . The exact solution (80)-(82) represents a one-parameter family of periodic waves of amplitude , traveling with the velocity v in a positive x-direction. The spatial period and time period for each wave are defined by
where is a complete elliptic integral of the first kind, see . The limit implies that the periodic behavior reduces to the propagation of a single solitary wave.
In order to compute the periodic case on approximately the same space-time domain as in the single solitary case, we again set the parameter values , , and the parameter k is experimentally set up as to have periods and .
The run of the algorithm is carried out up to one time period over the problem domain . The initial and nonhomogeneous Dirichlet conditions are extracted from the exact solution (80) and the same linear algebraic solver is used as in the previous case. Figure 2 depicts the propagation of approximation solutions of periodic waves from the initial condition to the final time for a piecewise linear approximation on the finest considered space-time mesh with time step and mesh size . Other coarse grids also produce similar plots.
Figure 2. The 3D plot of approximation solutions of periodic waves (left) and corresponding isolines in space-time domain (right).
In what follows, we shall proceed similarly as in Section 5.1 to verify the convergence and preservation of studied invariant quantities.
5.2.1 Convergence with respect to h
The h-convergence in the periodic case is investigated on a sequence of five successive refined grids partitioning the considered problem domain . The choice of time step is again small enough to suppress the influence of time discretization errors, and the computations are performed by piecewise linear approximations, subsequently.
The obtained results recorded in Tables 5 and 6 illustrate the same behavior of computational errors in the -norm and the energy norm with respect to the spatial discretization as in the case of a single solitary wave propagation. The computed EOCs at all four monitoring time instances keep asymptotically the same orders, i.e., and , for piecewise linear approximations and confirm the spatially suboptimal a priori error estimates (66) with respect to the -norm and spatially optimal estimates (67) in the energy norm, respectively. The influence of discretization errors on computations with a long time domain can be better eliminated by using the Crank-Nicolson numerical scheme instead of the backward Euler method.
Table 5. Periodic case: Computational errors in the -norm and experimental orders of convergence for approximation on a consequence of meshes at time instancest( )
Table 6. Periodic case: Computational errors in the energy norm and experimental orders of convergence for approximation on a consequence of meshes at time instancest( )
5.2.2 Convergence with respect to τ
The τ-convergence is experimentally verified by the computations on the finest spatial grid having 1,760 elements with piecewise linear approximation. The computations are performed by five different time steps τ and monitored at final time of one period . The theoretical results are in accordance with the observations listed in Table 7, i.e., and .
Table 7. Periodic case: Computational errors in the -norm and the energy norm for approximation with respect to time step ( )
From the presented numerical results in Sections 5.1.1-5.1.2 and 5.2.1-5.2.2, we see that the quality of approximate solutions obtained for a single solitary case and a periodic case is quite comparable.
5.2.3 Invariant conservation quantities
Similarly as in Section 5.2.3, we monitor the preservation of invariants of mass, momentum and energy defined by (75), (76) and (77), respectively. During the whole period of time, in the course of which the waves propagate inside the periodic domain , all these three invariants of motion remain conserved and equal to their original values that are well-determined analytically at .
The lack of similar problems in the literature caused that our experiments with periodic waves could not be compared with other methods, thus Table 8 captures only the development of errors in the -norm and the -norm and keeping the invariant quantities during the whole computation performed on the finest space-time grid. All three invariants of motion are not different from their analytical values, according to which this method can be considered suitable also for nonperiodic cases.
Table 8. Periodic case: Computed invariant quantities and errors in the -norm and the -norm ( , )
We have presented and theoretically analyzed an efficient numerical method for the solution of the RLW equation, which is based on the space dicretization by the discontinuous Galerkin method and a semi-implicit time discretization with suitable linearization of convective terms. Under some additional assumptions, we have derived a priori error estimates, namely in the -norm and in the energy norm. On the other hand, the presented numerical experiments for single solitary as well as periodic cases signal a better behavior of the experimental -order of convergence, which is expected to be asymptotically for piecewise linear approximations with a nonsymmetric variant of interior penalty Galerkin discretizations. In the case of the energy norm, we obtain the optimal experimental order of convergence.
The obtained results confirm that the proposed scheme is a powerful and reliable method for the numerical solution of a nonstationary nonlinear partial differential equation such as the RLW equation.
The authors declare that they have no competing interests.
Both authors contributed equally, read and approved the final version of the manuscript.
Dedicated to Professor Hari M Srivastava.
The authors would like to express their sincere gratitude to the referees for valuable comments and helpful suggestions. JH also would like to thank P. Červenková for her assistance with elaboration of numerical experiments. This work was partly supported by the ESF Project No. CZ.1.07/2.3.00/09.0155 ‘Constitution and improvement of a team for demanding technical computations on parallel computers at TU Liberec’ and by SGS Project ‘Modern numerical methods’ financed by TU Liberec.
Peregrine, DH: Calculations of the development of an undular bore. J. Fluid Mech.. 25, 321–330 (1966). Publisher Full Text
Haq, F, Islam, S, Tirmizi, IA: A numerical technique for solution of the MRLW equation using quartic B-splines. Appl. Math. Model.. 34(12), 4151–4160 (2010). Publisher Full Text
Islam, S, Haq, S, Ali, A: A meshfree method for the numerical solution of the RLW equation. J. Comput. Appl. Math.. 223(2), 997–1012 (2009). Publisher Full Text
Chen, Y, Mei, L: Explicit multistep method for the numerical solution of RLW equation. Appl. Math. Comput.. 218, 9547–9554 (2012). Publisher Full Text
Dag, I, Ozer, MN: Approximation of the RLW equation by the least square cubic B-spline finite element method. Appl. Math. Model.. 25, 221–231 (2001). Publisher Full Text
Esen, A, Kutluay, S: Application of a lumped Galerkin method to the regularized long wave equation. Appl. Math. Comput.. 174(2), 833–845 (2006). Publisher Full Text
Arnold, DN, Brezzi, F, Cockburn, B, Marini, LD: Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal.. 39(5), 1749–1779 (2002). Publisher Full Text
Dolejší, V, Feistauer, M, Hozman, J: Analysis of semi-implicit DGFEM for nonlinear convection-diffusion problems. Comput. Methods Appl. Mech. Eng.. 196, 2813–2827 (2007). Publisher Full Text
Dolejší, V, Feistauer, M, Sobotíková, V: Analysis of the discontinuous Galerkin method for nonlinear convection-diffusion problems. Comput. Methods Appl. Mech. Eng.. 194, 2709–2733 (2005). Publisher Full Text
Dolejší, V, Hozman, J: A priori error estimates for DGFEM applied to nonstationary nonlinear convection-diffusion equation. In: Kreiss G (ed.) Numerical Mathematics and Advanced Applications, pp. 459–468. Springer, Berlin ENUMATH, 2009. (2010)
Djidjeli, K, Price, WG, Twizell, EH, Cao, Q: A linearized implicit pseudo-spectral method for some model equations: the regularized long wave equations. Commun. Numer. Methods Eng.. 19, 847–863 (2003). Publisher Full Text