Research

# Analysis and application of the discontinuous Galerkin method to the RLW equation

Jiří Hozman1* and Jan Lamač12

Author Affiliations

1 Faculty of Science, Humanities and Education, Technical University of Liberec, Studentská 2, Liberec, 461 17, Czech Republic

2 Faculty of Mathematics and Physics, Charles University in Prague, Sokolovská 83, Prague, 186 75, Czech Republic

For all author emails, please log on.

Boundary Value Problems 2013, 2013:116  doi:10.1186/1687-2770-2013-116

 Received: 14 December 2012 Accepted: 23 April 2013 Published: 7 May 2013

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### Abstract

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

### 1 Introduction

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 [1]) 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 [2]. 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 [3], over collocation methods [4,5], to finite element approaches [6,7], or Galerkin methods [8], 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 [11]. Among several variants of DG methods, we deal with the nonsymmetric variant interior penalty Galerkin discretizations; see [12]. 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 [13]. Consequently, we extend the results from [13], and the attention is paid to the a priori error analysis of the method with the aid of standard techniques introduced in [14] and [15].

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 Ω = ( a , b ) R be a bounded open interval and T > 0 be a final time. We set Q T = Ω × ( 0 , T ) and consider the following initial boundary value problem: Find u ( x , t ) : Q T = Ω × ( 0 , T ) R such that, for all T > 0 ,

u t + u x + ε u u x μ t ( 2 u x 2 ) = 0 in  Q T , (1)

u ( a , t ) = u D a ( t ) and u ( b , t ) = u D b ( t ) for all  t ( 0 , T ) , (2)

u ( x , 0 ) = u 0 ( x ) for all  x Ω , (3)

where constant parameters ε > 0 and μ > 0 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 u 0 : Ω R . These given data have to satisfy the compatibility conditions prescribed at both endpoints of the domain Ω, i.e.,

u D a ( 0 ) = u 0 ( a ) and u D b ( 0 ) = u 0 ( b ) . (4)

The whole system (1)-(4) was found to have single solitary or periodic traveling wave solutions; for details, see [2].

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 k 0 be an integer and p [ 1 , ] . We use the well-known Lebesgue and Sobolev spaces L p ( Ω ) , H k ( Ω ) , Bochner spaces L p ( 0 , T ; X ) of functions defined in ( 0 , T ) with values in a Banach space X and the spaces C k ( [ 0 , T ] ; X ) of k-times continuously differentiable mappings of the interval [0, T] with values in X. By H 0 1 ( Ω ) we denote the subspace of all functions v H 1 ( Ω ) satisfying v ( a ) = v ( b ) = 0 . To this end, we use the following notation for a scalar product in L 2 ( Ω ) by

( u , v ) = Ω u v d x , u , v L 2 ( Ω ) (5)

for a norm in L 2 ( Ω ) by 2 = ( , ) , for a seminorm in H 1 ( Ω ) by | | 1 , 2 = ( , ) and for a norm in H 1 ( Ω ) by 1 , 2 = ( , ) + ( , ) . It is a known fact that | | 1 , 2 is a norm on H 0 1 ( Ω ) equivalent to 1 , 2

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 a ( , ) and a nonlinear convection form b ε ( , ) from [13], i.e.,

a ( u ( t ) , v ) = Ω u ( t ) x v d x , (6)

b ε ( u ( t ) , v ) = Ω f ( u ( t ) ) x v d x with  f ( u ) = u + ε 2 u 2 , (7)

where symbol u ( t ) stands for the function on Ω such that u ( t ) ( x ) , x Ω and function f ( u ) in (7) represents the physical flux.

Definition 1 We say that u is a weak solution of problem (1)-(3) if u L 2 ( 0 , T ; H 1 ( Ω ) ) L ( Q T ) such that u t L 2 ( 0 , T ; H 1 ( Ω ) ) and the following conditions are satisfied:

(a) u u L 2 ( 0 , T ; H 0 1 ( Ω ) ) , where  u ( t ) H 1 ( Ω )  such that (a) u ( t ) | x = a = u D a ( t )  and  u ( t ) | x = b = u D b ( t ) for a.e.  t ( 0 , T ) , (b) d d t ( u ( t ) , v ) + b ε ( u ( t ) , v ) + μ d d t a ( u ( t ) , v ) = 0 (b) v H 0 1 ( Ω )  and a.e.  t ( 0 , T ) , (c) ( u ( 0 ) , v ) = ( u 0 , v ) v H 0 1 ( Ω ) , u 0 L 2 ( Ω ) . (8)

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 H λ k ( ( a , b ) ) with period λ > 0 satisfying mod ( b a , λ ) = 0 .

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

Assumptions (R)

(9)
(10)
(11)

### 3 Discretization

Let T h ( h > 0 ) be a family of the partitions of the closure Ω ¯ = [ a , b ] of the domain Ω into N closed mutually disjoint subintervals I k = [ x k 1 , x k ] with length h k : = x k x k 1 and the symbol J stands for an index set { 1 , , N } . Then we call T h = { I k , k J } a partition with a spatial step h : = max k J ( h k ) and interval I k an element. By E h we denote the set of all partition nodes of Ω, i.e., E h = { x 0 = a , x 1 , , x N 1 , x N = b } . Further, we label by E h I the set of all inner nodes. Obviously, E h = E h I { a , b } .

We additionally assume that the partitions satisfy the following condition.

Assumption (M) T h are locally quasi-uniform:

C q 1 : h k C q h k I k , I k T h   sharing a node . (12)

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 s k N and a local polynomial degree p k N to each I k T h . Then we set the vectors s { s K , K T h } and p { p K , K T h } . Over the triangulation T h , we define the so-called broken Sobolev space corresponding to the vector s as

H s ( Ω , T h ) : = { v ; v | I k H s k ( I k ) I k T h } (13)

with the norm

v H s ( Ω , T h ) : = ( I k T h v H s k ( I k ) 2 ) 1 / 2 (14)

and the seminorm

| v | H s ( Ω , T h ) : = ( I k T h | v | H s k ( I k ) 2 ) 1 / 2 , (15)

where H s k ( I k ) and | | H s k ( I k ) denote the standard norm and the seminorm on the Sobolev space H s k ( I k ) , I k T h .

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

S h p S h p ( Ω , T h ) : = { v ; v | I k P p k ( I k ) I k T h } , (16)

where P p k ( I k ) denotes the space of all polynomials of degree p k on I k , I k T h .

Let us denote v ( x k + ) = lim ε 0 + v ( x k + ε ) and v ( x k ) = lim ε 0 + v ( x k ε ) . Then we can define the jump and average of v at inner points x k E h I of the domain Ω by

[ v ( x k ) ] = v ( x k ) v ( x k + ) , v ( x k ) = 1 2 ( v ( x k ) + v ( x k + ) ) . (17)

By convention, we also extend the definition of jump and mean value for endpoints of Ω, i.e., [ v ( x 0 ) ] = v ( x 0 + ) , v ( x 0 ) = v ( x 0 + ) , [ v ( x N ) ] = v ( x N ) and v ( x N ) = v ( x N ) . In case that x k E h are arguments of v ( x k ) or v ( x k + ) , we usually omit these arguments x k , x k + and write simply v and v + , respectively.

#### 3.1 DG semi-discrete formulation

Now, we recall the space semi-discrete DG scheme presented in [11]. First, we multiply (1) by a test function v h S h p , integrate over an element I k T h and use integration by parts in the dispersion term a h and convection term b h ε of (1) subsequently. Further, we sum over all I k T h and add some artificial terms vanishing for the exact solution such as penalty J h σ 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 [16]) for the discretization of the convection term and end up with the following DG formulation for the semi-discrete solution u h ( t ) , introduced in [13] as a system of ordinary differential equations, i.e.,

d d t { ( u h ( t ) , v h ) + μ a h ( u h ( t ) , v h ) + μ J h σ ( u h ( t ) , v h ) } + b h ε ( u h ( t ) , v h ) = 0 v h S h , t ( 0 , T ) , (18)

where forms a h ( , ) and b h ε ( , ) stand for the semi-discrete variants of forms (6) and (7), i.e.,

a h ( u ( t ) , v ) = k J I k u ( t ) x v d x x E h u ( t ) x [ v ] + x E h v [ u ( t ) ] , (19)

b h ε ( u ( t ) , v ) = k J I k ( u + ε 2 u 2 ) v d x + x E h H ( u ( t ) , u + ( t ) ) [ v ] . (20)

The crucial item of the DG formulation of the model problem is the treatment of the convection part. We proceed analogously as in [13], where the convection terms are approximated with the aid of the following numerical flux H ( , ) through node x E h in the positive direction (i.e., outer normal is equal to one):

H ( u ( x ) , u ( x + ) ) = { f ( u ( x ) ) = u ( x ) + ε 2 u 2 ( x ) , if  A > 0 , f ( u ( x + ) ) = u ( x + ) + ε 2 u 2 ( x + ) , if  A 0 , (21)

where A = f ( u ( x ) + u ( x + ) 2 ) and the choice of f ( u ( x 0 ) ) and f ( u ( x N + ) ) for boundary points has to satisfy the prescribed Dirichlet boundary conditions; for more details, see [14].

In what follows, we shall assume that the numerical flux H : R 2 R has the following properties.

Assumptions (H)

(H1) H ( u , v ) is Lipschitz-continuous with respect tou, v:

| H ( u , v ) H ( u , v ) | C H ( | u u | + | v v | ) , u , v , u , v R . (22)

(H2) H ( u , v ) is consistent:

H ( u , u ) = f ( u ) , u R . (23)

(H3) H ( u , v ) is conservative:

H ( u , v ) = H ( v , u ) in the negative direction , u , v R . (24)

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 R 2 .

A particular attention should be also paid to the treatment of the dispersion terms, which include an artificially added stabilization in the form x E h v [ u ( t ) ] , 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

J h σ ( u ( t ) , v ) = x E h I σ [ u ( t ) ] [ v ] + σ ( x 0 ) ( u ( x 0 + , t ) u D a ( t ) ) v ( x 0 + ) + σ ( x N ) ( u ( x N , t ) u D b ( t ) ) v ( x N ) (25)

which replaces the inter-element discontinuities and guarantees the fulfillment of the prescribed boundary conditions.

The penalty parameter function σ : E h R for a nonsymmetric variant is defined in spirit of [17] as

σ ( x ) = { p 1 2 / h 1 , x = a , min ( p k 2 / h k , p k + 1 2 / h k + 1 ) , x E h I { x } = I k I k + 1 , p N 2 / h N , x = b . (26)

In order to simplify the notation, we introduce the form

A h μ ( u ( t ) , v ) : = ( u ( t ) , v ) + μ a h ( u ( t ) , v ) + μ J h σ ( u ( t ) , v ) , u ( t ) , v S h p , t ( 0 , T ) , (27)

which is bilinear due to (19) and (25). Consequently, we can here define the semi-discrete solution u h of problem (8).

Definition 2 We say that u h is a semidiscrete solution of problem (8) if u h C 1 ( 0 , T ; S h p ) and the following conditions are satisfied:

(a) d d t A h μ ( u h ( t ) , v ) + b h ε ( u h ( t ) , v h ) = 0 v h S h p , t ( 0 , T ) , (b) ( u h ( 0 ) , v h ) = ( u 0 , v h ) v h S h p . (28)

#### 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 [13], we have proposed a semi-implicit time discretization based on the backward Euler scheme with the linearized convection form b h ε 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 [ 0 , T ] as 0 = t 0 < t 1 < t 2 < < t N = T , denoting each time step by τ l t l + 1 t l and letting u h l stand for the approximate solution of u h ( t l ) , t l [ 0 , T ] , l = 0 , , M . The linearization of the physical flux f is treated in spirit of [13] as

f ( u h l + 1 ) ( 1 + ε u h l ) u h l + 1 ε 2 ( u h l ) 2 , l = 0 , , M , (29)

which implies the splitting of a convection form in the following way:

b h ε ( u h l + 1 , v h ) b h L ε ( u h l , u h l + 1 , v h ) b h N ε ( u h l , v h ) (30)

with

b h L ε ( u h l , u h l + 1 , v h ) = k J I k ( 1 + ε u h l ) u h l + 1 v h d x + x E h H L ( ( u h l ) , ( u h l ) + , ( u h l + 1 ) , ( u h l + 1 ) + ) [ v h ] , (31)

b h N ε ( u h l , v h ) = k J I k ε 2 ( u h l ) 2 v h d x + x E h H N ( ( u h l ) , ( u h l ) + ) [ v h ] , (32)

where H L ( , , , ) and H N ( , ) represent the corresponding linearized and nonlinear parts of the original numerical flux H ( , ) given by (21); for more details, see [13]. One can easily observe that the form b h L ε ( , , ) is linear in its second and third argument and the form b h N ε ( , ) 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 0 = t 0 < t 1 < < t r = T be a partition of the time interval [ 0 , T ] and τ l t l + 1 t l , l = 0 , 1 , , M . We define the approximate solution of problem (8) as functions u h l S h p , t [ 0 , T ] , l = 0 , , M , satisfying the following conditions:

(a) A h μ ( u h l + 1 , v h ) + τ l b h L ε ( u h l , u h l + 1 , v h ) = A h μ ( u h l , v h ) + τ l b h N ε ( u h l , v h ) v h S h p , (b) u h 0  is  S h p  approximation of  u 0 . (33)

Discrete problem (33) is equivalent to a system of linear algebraic equations at each time instant t l [ 0 , T ] . 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 u h l S h p , τ l and t l + 1 [ 0 , T ] , we find u h l + 1 such that

A h l ( u h l + 1 , v h ) = f h l ( v h ) v h S h p , (34)

where

A h l ( u h l + 1 , v h ) = A h μ ( u h l + 1 , v h ) + τ l b h L ε ( u h l , u h l + 1 , v h ) , u h l , v h S h p (35)

f h l ( v h ) = A h μ ( u h l , v h ) + τ l b h N ε ( u h l , v h ) , v h S h p . (36)

Using the definitions (27) and (31), one can see that A h l is a bilinear form on the finite dimensional space and f h l is a linear functional. Moreover, the form A h l is coercive, i.e.,

A h l ( v h , v h ) v h 2 2 v h S h p . (37)

Hence, equation (34) has a unique solution u h l + 1 S h p . □

### 4 A priori error analysis

For error analysis and in experiments, we consider p k = p for all k J . Thus we denote S h p = S h p . Now we would like to analyze the error estimates of the approximate solution u h l , l = 0 , 1 ,  , obtained by method (33). For simplicity, we consider a uniform partition t l = l τ , l = 0 , 1 , , M , of the time interval [ 0 , T ] with time step τ = T / M , where M > 1 is an integer.

Let Π h u l be the standard S h p -interpolation of u l = u ( t l ) , ( l = 0 , , M ) satisfying (cf.[18]) for all v H p + 1 ( I k ) , I k T h ,

Π h v v L 2 ( I k ) c ˜ h k p + 1 | v | H p + 1 ( I k ) , (38)

| Π h v v | H 1 ( I k ) c ˜ h k p | v | H p + 1 ( I k ) , (39)

for a generic constant c ˜ > 0 independent of h and v. We set

ξ h l = u h l Π h u l S h p , η h l = Π h u l u l H p + 1 ( Ω , T h ) . (40)

Then the error e h l = u h l u l can be expressed as

e h l = ξ h l + η h l , l = 0 , 1 , , M . (41)

Setting ξ h l + 1 in (33), we get

A h μ ( u h l + 1 u h l , ξ h l + 1 ) + τ ( b h L ε ( u h l , u h l + 1 , ξ h l + 1 ) b h N ε ( u h l , ξ h l + 1 ) ) = 0 . (42)

On the other hand, from (18) it follows

A h μ ( u t ( t l + 1 ) , ξ h l + 1 ) + b h ε ( u ( t l + 1 ) , ξ h l + 1 ) = 0 . (43)

Multiplying (43) by τ, subtracting from (42) and using again the linearity of the form A h μ , we get

A h μ ( u h l + 1 u h l τ u t ( t l + 1 ) , ξ h l + 1 ) + τ ( b h L ε ( u h l , u h l + 1 , ξ h l + 1 ) b h N ε ( u h l , ξ h l + 1 ) b h ε ( u ( t l + 1 ) , ξ h l + 1 ) ) = 0 . (44)

Since u h l + 1 u h l = ξ h l + 1 + η h l + 1 + u l + 1 u l ξ h l η h l , we can rewrite equation (44) in the following way:

A h μ ( ξ h l + 1 ξ h l , ξ h l + 1 ) = A h μ ( u l + 1 u l τ u t ( t l + 1 ) , ξ h l + 1 ) A h μ ( η h l + 1 η h l , ξ h l + 1 ) + τ ( b h L ε ( u h l , u h l + 1 , ξ h l + 1 ) b h N ε ( u h l , ξ h l + 1 ) b h ε ( u ( t l + 1 ) , ξ h l + 1 ) ) . (45)

For the term on the right-hand side of equation (45), we use decomposition ( m n ) m = 1 2 ( m 2 n 2 + ( m n ) 2 ) and linearity of the form A h μ . Together with (30), we finally get

1 2 { A h μ ( ξ h l + 1 , ξ h l + 1 ) A h μ ( ξ h l , ξ h l ) + A h μ ( ξ h l + 1 ξ h l , ξ h l + 1 ξ h l ) } = A h μ ( u l + 1 u l τ u t ( t l + 1 ) , ξ h l + 1 ) A h μ ( η h l + 1 η h l , ξ h l + 1 ) + τ { b h L ε ( u h l , u h l + 1 , ξ h l + 1 ) b h L ε ( u h l , u ( t l + 1 ) , ξ h l + 1 ) + b h L ε ( u h l , u ( t l + 1 ) , ξ h l + 1 ) b h L ε ( u ( t l + 1 ) , u ( t l + 1 ) , ξ h l + 1 ) + b h N ε ( u ( t l + 1 ) , ξ h l + 1 ) b h N ε ( u h l , ξ h l + 1 ) } . (46)

For next estimates, we use the following lemmas.

Lemma 2Under assumptions (R) for t l , t l + 1 [ 0 , T ] , the following hold:

| ( u l + 1 u l τ u t ( t l + 1 ) , ξ h l + 1 ) | c τ 2 ξ h l + 1 2 , (47)

| ( η h l + 1 η h l , ξ h l + 1 ) | c τ h 2 ( p + 1 ) ξ h l + 1 2 , (48)

| a h ( u l + 1 u l τ u t ( t l + 1 ) , ξ h l + 1 ) | c τ 2 ( | ξ h l + 1 | 1 , 2 + J h σ ( ξ h l + 1 , ξ h l + 1 ) 1 / 2 ) , (49)

| a h ( η h l + 1 η h l , ξ h l + 1 ) | c τ h 2 ( p + 1 ) ( | ξ h l + 1 | 1 , 2 + J h σ ( ξ h l + 1 , ξ h l + 1 ) 1 / 2 ) , (50)

| J h σ ( u l + 1 u l τ u t ( t l + 1 ) , ξ h l + 1 ) | c τ 2 J h σ ( ξ h l + 1 , ξ h l + 1 ) 1 / 2 , (51)

| J h σ ( η h l + 1 η h l , ξ h l + 1 ) | c τ h 2 ( p + 1 ) J h σ ( ξ h l + 1 , ξ h l + 1 ) 1 / 2 , (52)

wherecis a generic constant independent ofhandτ.

Proof The proof of these standard estimates can be found, for instance, in [15]. □

Lemma 3Under assumptions (R), (H) and for t l , t l + 1 [ 0 , T ] , the following hold:

| b h L ε ( u h l , u h l + 1 u ( t l + 1 ) , ξ h l + 1 ) | c ( | ξ h l + 1 | 1 , 2 + J h σ ( ξ h l + 1 , ξ h l + 1 ) 1 / 2 ) ( h p + 1 + ξ h l + 1 2 ) , (53)

| b h L ε ( u h l , u ( t l + 1 ) , ξ h l + 1 ) b h L ε ( u ( t l + 1 ) , u ( t l + 1 ) , ξ h l + 1 ) | c ( | ξ h l + 1 | 1 , 2 + J h σ ( ξ h l + 1 , ξ h l + 1 ) 1 / 2 ) ( h p + 1 + ξ h l 2 + τ ) , (54)

| b h N ε ( u ( t l + 1 ) , ξ h l + 1 ) b h N ε ( u h l , ξ h l + 1 ) | c ( | ξ h l + 1 | 1 , 2 + J h σ ( ξ h l + 1 , ξ h l + 1 ) 1 / 2 ) ( h p + 1 + ξ h l 2 + τ ) , (55)

wherecis a generic constant independent ofhandτ.

Proof Again, one can find the proof of these estimates in [15]. □

Since A h μ ( ξ h l + 1 ξ h l , ξ h l + 1 ξ h l ) is always nonnegative, applying previous lemmas gives us

1 2 A h μ ( ξ h l + 1 , ξ h l + 1 ) 1 2 A h μ ( ξ h l , ξ h l ) c { τ ( τ + h p + 1 ) ξ h l + 1 2 + μ τ ( τ + h p ) ( | ξ h l + 1 | 1 , 2 + J h σ ( ξ h l + 1 , ξ h l + 1 ) 1 / 2 ) + μ τ ( τ + h p ) J h σ ( ξ h l + 1 , ξ h l + 1 ) 1 / 2 + τ ( | ξ h l + 1 | 1 , 2 + J h σ ( ξ h l + 1 , ξ h l + 1 ) 1 / 2 ) ( ξ h l + 1 2 + h p + 1 + ξ h l 2 + τ ) } . (56)

Multiplying by 2, applying the Young inequality and using the definition of the form A h μ , we obtain

ξ h l + 1 2 2 + μ ( | ξ h l + 1 | 1 , 2 2 + J h σ ( ξ h l + 1 , ξ h l + 1 ) ) ξ h l 2 2 + μ ( | ξ h l | 1 , 2 2 + J h σ ( ξ h l , ξ h l ) ) + 2 c τ { τ 2 2 ν 1 + ν 1 2 ξ h l + 1 2 2 + h 2 ( p + 1 ) 2 ν 2 + ν 2 2 ξ h l + 1 2 2 + μ ( τ 2 ν 3 + h 2 ( p ) ν 4 + ν 3 + ν 4 2 ( | ξ h l + 1 | 1 , 2 2 + J h σ ( ξ h l + 1 , ξ h l + 1 ) ) ) + μ ( τ 2 2 ν 5 + h 2 ( p ) 2 ν 6 + ν 5 + ν 6 2 J h σ ( ξ h l + 1 , ξ h l + 1 ) ) + ν 7 + ν 8 + ν 9 + ν 10 2 μ ( | ξ h l + 1 | 1 , 2 2 + J h σ ( ξ h l + 1 , ξ h l + 1 ) ) + 1 μ ν 7 ξ h l + 1 2 2 + h 2 ( p + 1 ) μ ν 8 + 1 μ ν 9 ξ h l 2 2 + τ 2 μ ν 10 } . (57)

If we take into account that J h σ ( ξ h l + 1 , ξ h l + 1 ) | ξ h l + 1 | 1 , 2 2 + J h σ ( ξ h l + 1 , ξ h l + 1 ) , the previous inequality can be rewritten as

( 1 c τ ( ν 12 + 2 μ ν 7 ) ) ξ h l + 1 2 2 + ( 1 c τ i = 3 10 ν i ) μ ( | ξ h l + 1 | 1 , 2 2 + J h σ ( ξ h l + 1 , ξ h l + 1 ) ) ( 1 + 2 c τ ν 9 μ ) ξ h l 2 2 + μ ( | ξ h l | 1 , 2 2 + J h σ ( ξ h l , ξ h l ) ) + c τ q ( τ , h , μ ) , (58)

where we denoted ν 12 = ν 1 + ν 2 and

q ( τ , h , μ ) = ( 1 ν 1 + 2 μ ν 3 + μ ν 5 + 2 μ ν 10 ) τ 2 + ( 2 ν 4 + 1 2 ν 6 ) μ h 2 p + ( 1 ν 2 + 2 ν 8 μ ) h 2 ( p + 1 ) . (59)

Let us now introduce the so-called energy norm

v h = | v h | 1 , 2 2 + J h σ ( v h , v h ) (60)

and the norm

v h μ = v h 2 2 + μ v h 2 . (61)

Denoting C L = c max { ν 12 + 2 μ ν 7 , i = 3 10 ν i } and C R = 2 c τ μ ν 9 from (58), it follows

( 1 τ C L ) ξ h l + 1 μ 2 ( 1 + C R ) ξ h l μ 2 + c τ q ( τ , h , μ ) . (62)

In order to finish our estimates, we require a fulfillment of the following technical assumption.

Assumption (T)

(T1) There exists θ ( 0 , 1 ) such that 0 < τ < θ / C L .

If assumption (T) is fulfilled, then τ < 1 C L μ μ ν 12 + 2 / ν 7 ν 7 2 μ . Thus we can also reformulate assumption (T) so that τ = O ( μ ) .

Thus, let us assume that assumption (T) holds, then

ξ h l + 1 μ 2 B ξ h l μ 2 + c τ 1 τ C L q ( τ , h , μ ) (63)

with B = 1 + τ C R 1 τ C L = 1 + τ C L + C R 1 τ C L exp ( τ C L + C R 1 τ C L ) . Consequently,

ξ h l μ 2 B l ξ h 0 μ 2 + B l 1 B 1 c τ 1 τ C L q ( τ , h , μ ) (64)

and since B 1 = τ C L + C R 1 τ C L , we have c τ ( B 1 ) ( 1 τ C L ) = c C L + C R , i.e.,

ξ h l μ 2 B l ξ h 0 μ 2 + ( B l 1 ) c C L + C R q ( τ , h , μ ) exp ( l τ C L + C R 1 τ C L ) ( ξ h 0 μ 2 + c C L + C R q ( τ , h , μ ) ) C ˜ exp ( T C L + C R 1 τ C L ) ( μ h 2 p + h 2 ( p + 1 ) + μ ν 7 ν 9 2 ( ν 7 + ν 9 ) q ( τ , h , μ ) ) , (65)

where we used a straightforward estimate ξ h 0 μ 2 C ˜ ( μ h 2 p + h 2 ( p + 1 ) ) . We can notice that due to the presence of the factor μ in front of the function q ( τ , h , μ ) on the left-hand side of (65), we lost μ in denominators of q ( τ , h , μ ) .

Now we are ready to formulate the main theorem.

Theorem 1Let assumptions (M), (H), (R) and (T) be satisfied, then there exists a constant C = C ( μ ) such that

(66)

(67)

where is defined by (60).

Proof Since e h l μ ξ h l μ + η h l μ , the statement of the theorem comes from (65) and the fact that η h l μ c ˜ ( h p + 1 + μ h p ) . Then we set

C ( μ ) = C ˜ exp ( T C L + C R 1 θ ) ν 7 ν 9 ν 7 + ν 9 max { 1 ν j , j { 1 , 2 , , 6 , 8 , 10 } } . (68)

□

Remark 3 Theorem 1 implies that the error of our method is O ( h p + τ ) in both energy and L 2 -norm. However, as we will see in the next section, the error estimate in the L 2 -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 [19]. This is useful namely in the cases when convection terms dominate, i.e., μ 0 + . Consequently, in these cases assumption (T) can be weakened to a CFL-like condition τ = O ( h α ) for suitable α > 0 .

### 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 L 2 -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 t = l τ during all computations in the corresponding norms, i.e.,

err h , τ 0 , l u h l u ( l τ ) 2 , (69)

err h , τ 1 , l u h l u ( l τ ) , (70)

where u ( l τ ) is a prescribed exact solution at time and u h l is the numerical solution at time level 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

err h , τ n = err h n + err τ n , n = 0 , 1 , (71)

where

err h n D n ˜ h a n , err τ n D n ˆ τ b n , n = 0 , 1 . (72)

The constants D n ˜ , n = 0 , 1 , are independent of τ and D n ˆ , n = 0 , 1 , are independent of h. The values a n , b n , n = 0 , 1 , are the orders of accuracy of the method in the corresponding considered norms. We define the experimental order of convergence (EOC) by

a n = log ( err h 1 , τ n / err h 2 , τ n ) log ( h 1 / h 2 ) and b n = log ( err h , τ 1 n / err h , τ 2 n ) log ( τ 1 / τ 2 ) , n = 0 , 1 . (73)

#### 5.1 Single solitary case

The RLW equation has the following analytical single solitary wave solution given by

u ( x , t ) = 3 c sech 2 ( B 0 ( x x ¯ v t ) ) with  B 0 = 1 2 ε c μ ( 1 + ε c ) , (74)

which represents a single solitary wave of amplitude 3c, traveling with the velocity v = 1 + ε c in a positive x-direction and located initially at the point x ¯ . 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 [5], we set the parameter values c = 0.1 , x ¯ = 0.0 , ε = μ = 1.0 . The run of the algorithm is carried out up to time T = 20.0 over the problem domain [ 40 , 60 ] . 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 τ = 0.001 and mesh size h = 0.05 . 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 τ = 10 3 . Numerical experiments are carried out with the use of piecewise linear ( P 1 ) 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 L 2 -norm and the energy norm at four time instances t = 5.0 , t = 10.0 , t = 15.0 and t = 20.0 , the corresponding EOC during all computations. Since the exact solution u ( t ) is sufficiently regular over Ω, it follows from Remark 3 that the theoretical error estimates are of order O ( h p + τ ) . 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 L 2 -norm, which is expected to be asymptotically O ( h 2 ) for piecewise linear ( p = 1 ) approximations. These observations also correspond with the finite element approach from [6], where the same example was studied.

Table 1. Single solitary case: Computational errors in the L 2 -norm and experimental orders of convergence for P 1 approximation on a consequence of meshes at time instancest( τ = 10 3 )

Table 2. Single solitary case: Computational errors in the energy norm and experimental orders of convergence for P 1 approximation on a consequence of meshes at time instancest( τ = 10 3 )

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 O ( h p ) . Finally, both estimates in Theorem 1 confirm the well-know attribute of DG schemes from the class of convection-diffusion problems, cf.[14] and [15].

#### 5.1.2 Convergence with respect to τ

Secondly, we verify experimentally the convergence of the method in the L 2 -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 t = T in the L 2 -norm and the energy norm, respectively. We observe that both computational errors have EOC of order O ( τ ) in the corresponding norms, which is again in a good agreement with derived theoretical results.

Table 3. Single solitary case: Computational errors in the L 2 -norm and the energy norm for P 1 approximation with respect to time step ( h = 0.05 )

#### 5.1.3 Invariant conservation quantities

Similarly as in [13], we shall monitor the three conservation quantities for the propagation of the single solitary wave corresponding to mass

I M ( u ) = Ω u d x , (75)

momentum

I P ( u ) = Ω ( u 2 + μ ( u ) 2 ) d x , (76)

and energy

I E ( u ) = Ω ( u 3 + 3 u 2 ) d x , (77)

with respect to the run of the proposed algorithm. The analytical values for the invariants on the entire real domain are given (in [20]) by

I M ( u ) = 6 c B 0 , I P ( u ) = 12 c 2 B 0 + 48 B 0 c 2 μ 5 , I E ( u ) = 36 c 2 B 0 + 144 c 3 5 B 0 . (78)

Moreover, for the purpose of a more accurate comparison with reference results, we introduce the discrete l -norm defined by

err , l max k J ( | u h l ( x k 1 + ) u ( x k 1 + , l τ ) | , | u h l ( x k ) u ( x k , l τ ) | ) (79)

assessing the accuracy of the method by measuring the difference between the numerical and analytic solutions u h and u, respectively.

Table 4 records the invariant quantities together with errors err , err 0 computed on the finest space-time grid and compares obtained results with several previously presented schemes given in [6,7,20] and [5]. 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 [5]. Moreover, we append the recent results from [13] 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 l -norm and the L 2 -norm ( h = 0.05 , τ = 10 3 )

#### 5.2 Periodic case

The family of periodic solutions of the RLW equation may be analytically written as (cf.[20])

u ( x , t ) = A 0 + A 1 dn 2 ( B 1 ( x x ¯ v t ) , k ) with  v = 1 + ε c (80)

and the parameters A 0 , A 1 and B 1 given by

A 0 = c ( 1 2 k 2 k 4 k 2 + 1 ) , A 1 = 3 c k 4 k 2 + 1 , (81)

B 1 = 1 2 ε c μ ( 1 + ε c ) 1 k 4 k 2 + 1 4 , (82)

where dn ( , k ) is the Jacobi elliptic function and k [ 0 , 1 ) stands for the elliptic modulus; for definitions and other properties, see [21]. The exact solution (80)-(82) represents a one-parameter family of periodic waves of amplitude A 0 + A 1 , traveling with the velocity v in a positive x-direction. The spatial period ω k and time period T k for each wave are defined by

ω k = 2 K ( k ) / B 1 and T k = ω k / v , (83)

where K ( k ) is a complete elliptic integral of the first kind, see [21]. The limit k 1 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 c = 0.1 , x ¯ = 0.0 , ε = μ = 1.0 and the parameter k is experimentally set up as k = 0.63048 to have periods T k 20.0 and ω k 22.0 .

The run of the algorithm is carried out up to one time period T k over the problem domain [ 2 ω k , 2 ω k ] . 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 T k for a piecewise linear approximation on the finest considered space-time mesh with time step τ = 0.001 and mesh size h = 0.05 . 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 [ 44 , 44 ] . 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 L 2 -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., err h 0 = O ( h 2 ) and err h 1 = O ( h ) , for piecewise linear approximations and confirm the spatially suboptimal a priori error estimates (66) with respect to the L 2 -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 L 2 -norm and experimental orders of convergence for P 1 approximation on a consequence of meshes at time instancest( τ = 10 3 )

Table 6. Periodic case: Computational errors in the energy norm and experimental orders of convergence for P 1 approximation on a consequence of meshes at time instancest( τ = 10 3 )

#### 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 T k . The theoretical results are in accordance with the observations listed in Table 7, i.e., err τ 0 = O ( h ) and err τ 1 = O ( h ) .

Table 7. Periodic case: Computational errors in the L 2 -norm and the energy norm for P 1 approximation with respect to time step ( h = 0.05 )

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 [ 2 ω k , 2 ω k ] , all these three invariants of motion remain conserved and equal to their original values that are well-determined analytically at t = 0 .

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 l -norm and the L 2 -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 l -norm and the L 2 -norm ( h = 0.05 , τ = 10 3 )

### 6 Conclusion

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 O ( h p + τ ) in the L 2 -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 ( L 2 ) -order of convergence, which is expected to be asymptotically O ( h 2 + τ ) 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.

### Competing interests

The authors declare that they have no competing interests.

### Authors’ contributions

Both authors contributed equally, read and approved the final version of the manuscript.

### Acknowledgements

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.

### References

1. Peregrine, DH: Calculations of the development of an undular bore. J. Fluid Mech.. 25, 321–330 (1966). Publisher Full Text

2. Pava, JA: Nonlinear Dispersive Equations: Existence and Stability of Solitary and Periodic Travelling Wave Solutions, Am. Math. Soc., Providence (2009)

3. Esen, A, Kutluay, S: A finite difference solution of the regularized longwave equation. Math. Probl. Eng.. 2006, (2006) Article ID 85743

4. 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

5. 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

6. Chen, Y, Mei, L: Explicit multistep method for the numerical solution of RLW equation. Appl. Math. Comput.. 218, 9547–9554 (2012). Publisher Full Text

7. 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

8. 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

9. Cockburn, B: Discontinuous Galerkin methods for convection dominated problems. In: Barth TJ, Deconinck H (eds.) High-Order Methods for Computational Physics, pp. 69–224. Springer, Berlin (1999)

10. Cockburn B, Karniadakis GE, Shu C-W (eds.): Discontinuous Galerkin Methods, Springer, Berlin (2000)

11. Rivière, B: Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation. Frontiers in Applied Mathematics, SIAM, Philadelphia (2008)

12. 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

13. Hozman, J: Discontinuous Galerkin method for numerical solution of the regularized long wave equation. AIP Conf. Proc.. 1497, 118–125 (2012)

14. 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

15. 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

16. Feistauer, M, Felcman, J, Straškraba, I: Mathematical and Computational Methods for Compressible Flow, Oxford University Press, Oxford (2003)

17. 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)

18. Ciarlet, PG: The Finite Elements Method for Elliptic Problems, North-Holland, Amsterdam (1979)

19. Kučera, V: On ε-uniform error estimates for singularly perturbed problems in the DG method. Numerical Mathematics and Advanced Applications 2011. 368–378 (2013)

20. 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

21. Abramowitz, M, Stegun, IA: Handbook of Mathematical Functions, Dover, New York (1965)