# Energy preserving moving mesh methods applied to the BBM equation

###### Abstract

Energy preserving numerical methods for a certain class of PDEs are derived, applying the partition of unity method. The methods are extended to also be applicable in combination with moving mesh methods by the rezoning approach. These energy preserving moving mesh methods are then applied to the Benjamin–Bona–Mahony equation, resulting in schemes that exactly preserve an approximation to one of the Hamiltonians of the system. Numerical experiments that demonstrate the advantages of the methods are presented.

## 1 Introduction

Numerical solutions of differential equations by standard methods will typically not inherit invariant properties from the original, continuous problem. Since the energy-preserving methods of Courant, Friedrichs and Lewy were introduced in [1], the development of conservative methods has garnered much interest and considerable research, surveyed in [2] up to the early 1990s. In some important cases, conservation properties can be used to ensure numerical stability or existence and uniqueness of the numerical solution. In other cases, the conservation of one or more invariants can be of importance in its own right. In addition, as noted in [3], one may expect that when properties of the continuous dynamical system are inherited by the discrete dynamical system, the numerical solution can be more accurate, especially over large time intervals.

The discrete gradient methods for ordinary differential equations (ODEs), usually attributed to Gonzalez [4], are methods that preserve first integrals exactly. Since the late 1990s, a number of researchers have worked on extending this theory to create a counterpart for partial differential equations (PDEs), see e.g. [5, 6]. Such methods, which are either called discrete variational derivative methods or discrete gradient methods for PDEs, aim at preserving some discrete approximation of a first integral which is preserved by the continuous system. Up to very recently, the schemes presented have typically been based on a finite difference approach, and exclusively on fixed, uniform grids. Two different discrete variational derivative methods on fixed, non-uniform grids were presented by Yaguchi, Matsuo and Sugihara in [7, 8]. In [9], Miyatake and Matsuo introduce integral preserving methods on adaptive grids for certain classes of PDEs. Eidnes, Owren and Ringholm presented in [10] a general approach to extending the theory of discrete variational derivative methods, or discrete gradient methods for PDEs, to adaptive grids, using either a finite difference approach, or the partition of unity method, which can be seen as a generalization of the finite element method.

In this paper, we present an application of the approach introduced in [10] to the Benjamin–Bona–Mahony (BBM) equation, also called the regularized long wave equation in the literature. Although what we present here is a finite element method, the theory can be easily applied in a finite difference setting. Previously, there have been developed integral preserving methods for this equation [11], as well as adaptive moving mesh methods [12], but the schemes we are to present here are, to our knowledge, the first combining these properties. In fact, in [12] it is noted that combining integral preservation with adaptivity is an interesting topic for further research.

## 2 THE DISCRETE GRADIENT METHODS FOR PDEs

We give a quick survey of the discrete gradient methods for PDEs, and present an approach to the spatial discretization by the partition of unity method (PUM).

### 2.1 Problem statement

Consider a PDE of the form

(1) |

where denotes itself and its partial derivatives of any order with respect to the spatial variables , and where we assume that is sufficiently regular to allow all operations used in the following.

We define a first integral of (1) to be a functional satisfying

recalling that the variational derivative is defined as the function satisfying

This means that is preserved over time by (1), since

Furthermore, we may observe that if there exists some operator , skew-symmetric with respect to the inner product, such that

then is a first integral of (1), and we can state (1) on the form

(2) |

The idea behind the discrete variational derivative methods is to derive a discrete version of the PDE on the form (2), by obtaining a so-called discrete variational derivative and approximate by a skew-symmetric matrix, see e.g. [5].

As proven in [10], all discrete variatonal derivative methods can be expressed as discrete gradient methods on a system of ODEs obtained by discretizing (2) in space, to get a system

(3) |

where is a skew-symmetric matrix. The discrete gradient methods for such a system of ODEs preserve the first integral [13]. These numerical methods are given by

where is a consistent skew-symmetric time-discrete approximation to and is a discrete gradient of , defined as a function satisfying

There are many possible choices of discrete gradients. For the numerical experiments in this note, we will use the Average Vector Field (AVF) discrete gradient [6], given by

Note that when discretizing the system (2) in space, we do so by finding a discrete approximation to the integral , and define an energy preserving method to be a method preserving this approximation.

### 2.2 Partition of unity method on a fixed mesh

The partition of unity method is a generalization of the finite element method (FEM). Stating a weak form of (2), the problem consists of finding such that

We define an approximation to by

where the test functions span a finite-dimensional subspace . Referring to [10] for details, we then obtain the Galerkin form of the problem: Find such that

where, with ,

We end up with an ODE for the coefficients :

(4) |

Here, is a skew-symmetric matrix, with given by , and the system is thereby of the form (3). Then, the scheme

will preserve the approximated first integral in the sense that .

## 3 Adaptive Schemes

The primary motivation for using an adaptive mesh is usually to increase accuracy while keeping computational cost low, by improving discretization locally. Such methods are typically useful for problems with e.g. traveling wave solutions and boundary layers. The different strategies for adaptive meshes can be classified into two main groups [14]: The quasi-Lagrange approach involves coupling the evolution of the mesh with the PDE, and then solving the problems simultaneously; The rezoning approach consists of calculating the function values and mesh points in an intermittent fashion. Our method can be coupled with any adaptive mesh strategy utilizing the latter approach.

### 3.1 Adaptive discrete gradient methods

Let , , , and denote the discretization parameters and the numerical values obtained at the current time step and next time step, respectively. Note that we now alter the notion of a preserved first integral further, to requiring that . The idea behind our approach is to find based on and , transfer to to obtain , and then use to propagate in time to get . If the transfer operation between the meshes is preserving, i.e. if , then the next time step can be taken with the discrete gradient method for static meshes. If, however, non-preserving transfer is used, corrections are needed in order to get a numerical scheme. We introduce in [10] the scheme

(5) |

where is a vector which should be chosen so as to obtain a minimal correction, and such that . In the numerical experiments to follow, we have used .

A preserving transfer can by obtained using the method of Lagrange multipliers. Depending on whether - - or -refinement (or a combination) is used between time steps, we expect the shape and/or number of basis functions to change. See e.g. [14] or [15] for examples of how the basis may change through adaptivity. Denote by the trial space from the current time step and by the trial space for the next time step, and note that in general, . We wish to transfer the approximation from to , obtaining an approximation , while conserving the first integral, i.e. . This can be formulated as a constrained minimization problem:

(6) |

Observe that

where , and . The problem (6) can thus be reformulated as

This is a quadratic minimization problem with one nonlinear equality constraint, for which the solution is the solution of the nonlinear system of equations

which can be solved numerically using a suitable nonlinear solver.

## 4 Adaptive Energy Preserving Schemes for the Bbm Equation

### 4.1 The BBM equation

The BBM equation was introduced by Peregrine [16], and later studied by Benjamin et al. [17] as a model for small amplitude long waves on the surface of water in a channel. Conservative finite difference schemes for the BBM equation were proposed in [18] and [11], the latter being a discrete gradient method on fixed grids. A moving mesh FEM scheme employing a quasi-Lagrange approach is presented by Lu, Huang and Qiu in [12], which we also refer to for a more extensive list of references to the existing numerical schemes for the BBM equation.

Consider now an initial-boundary value problem of the one-dimensional BBM equation with periodic boundary conditions,

(7) | ||||

(8) | ||||

(9) |

By introducing the new variable , equation (7) can be rewritten on the form (2) as

for two different pairs of an antisymmetric differential operator and a Hamiltonian :

and

### 4.2 Discrete schemes

We apply the PUM approach to create numerical schemes which preserve an approximation to either or , splitting into elements . Defining the matrices and by their components

we set Note that the matrices and depend on the mesh, and thus will change when adaptivity is used. We will then distinguish between matrices from different time steps by writing e.g. and .

Approximating by as in section 2.2, we find

The integrals can be evaluated exactly and efficiently by considering elementwise which basis functions are supported on the element before applying Gaussian quadrature to obtain exact evaluations of the polynomial integrals. We define the matrix by

An approximation to the gradient of with respect to is found by the AVF discrete gradient

Thus we have the required terms for forming the system (4) and applying the adaptive discrete gradient method to it. Corresponding to (5), we get the scheme

Here we have chosen the skew-symmetric matrix to be a function of and , but could also have chosen e.g. , resulting in a decreased computational cost at the expense of less precise results. During testing, the basis functions were chosen as piecewise cubic polynomials.

In the same manner we may obtain a scheme that preserves . In this case

and

Note that the skew-symmetric matrix is independent of .

Defining the tensor by its elements

we get, with the convention of summation over repeated indices, the AVF discrete gradient with respect to given by the elements

and again the discrete gradient with respect to by

If we employ integral preserving transfer between the meshes, we get the scheme

where we note that is a skew-symmetric matrix. If non-preserving transfer is used, we need a correction term, as in the scheme above. The calculation of such a term is straightforward, but we omit it here for reasons of brevity.

To approximate the third derivative in , we need basis functions of at least degree three, and to guarantee skew-symmetry in , these basis functions need to be on the element boundaries. This is not obtainable with regular nodal FEM basis functions, so we have instead used third order B-spline basis functions as described in [19] during testing.

## 5 Numerical Results

To demonstrate the performance of our methods, we have tested them on two one-dimensional simple problems: A soliton solution, and the interaction of two waves. We have tested our - and -preserving schemes on uniform and moving meshes, and compared the results to those obtained using the explicit midpoint method. For the transfer operation between meshes, we have used a piecewise cubic interpolation method in the preserving scheme, and exact transfer in the preserving scheme.

### 5.1 Mesh adaptivity

As noted in section 3, our methods can be coupled with any adaptive mesh strategy using the rezoning approach. For our numerical experiments, we have used a simple method for -adaptivity based on the equidistribution principle: Splitting into intervals, we require that

where the monitor function is a function measuring how densely grid points should lie, based on the value of . For a general discussion on the choice of an optimal monitor function, see e.g. [20, 21]. For the problems we have studied, a generalized solution arc length monitor function proved to yield good results. This is given by

For , this is the usual arc length monitor function, in which case the equidistribution principle amounts to requiring that the arc length of over each interval is equal. In applications, we only have an approximation of , and hence must be approximated as well. We have applied a finite difference approximation and obtained approximately equidistributing grids using de Boor’s method as explained in [14, pp. 36-38].

### 5.2 Soliton solution

With , the exact solution of (7)–(9) is

with . This is a soliton solution which travels with a constant speed in -direction while maintaining its initial shape.

To evaluate the numerical solutions, we have compared them to the exact solution and calculated errors in shape and phase. The phase error is evaluated as

where , i.e. the location of the peak of the soliton in the numerical solution. The shape error is given by

where the peak of the exact solution is translated to match the peak of the numerical solution, and the difference in the shapes of the solitons is calculated.

The results of the numerical tests can be seen in figures 1–3. Here, denotes the degrees of freedom used in the spatial approximation and the fixed time step size. DG1 and DG1MM denotes the preserving scheme with fixed, uniform grid and adaptive grid, respectively; similary DG2 and DG2MM denotes the preserving scheme with uniform and adaptive grids.

In Figure 1 we see the relative errors in and . The DG1 and DG1MM schemes are compared to schemes using the same 3rd order nodal basis functions, but the trapezoidal rule for time-stepping, denoted by TR and TRMM. Likewise, the DG2 and DG2MM schemes are compared to the IM and IMMM schemes, using B-spline basis functions and the implicit midpoint method for discretization in time. The error in is very small for the DG1 and DG1MM schemes, as expected. Also the error in is very small for the DG2 and DG2MM schemes. The order of the error is not machine precision, but is instead dictated by the precision with which the nonlinear equations in each time step is solved. We can also see that while the TR and IM schemes, with and without moving meshes, have poor conservation properties, the moving mesh DG schemes seem to preserve quite well even the integrals they are not designed to preserve.

In figures 2 and 3 we see the phase and shape errors, of our methods compared to non-moving mesh methods and non-preserving methods, respectively. The advantage of using moving meshes is clear, especially for the preserving schemes. The usefulness on integral preservation is ambiguous in this case. It seems that what we gain in precision in phase, we lose in precision in shape, and vice versa.

### 5.3 A small wave overtaken by a large one

A typical test problem for the BBM equation is the interaction between two solitary waves. With an initial condition

one wave will eventually be overtaken by the other as long as , i.e. if one wave is larger than the other. There is no available analytical solution for this problem. The two waves are not solitons, as the amplitudes will change a bit after the waves have interacted [22].

Solutions obtained by solving the problem with our two energy preserving schemes, giving very similar results, are plotted in Figure 4. Also, to illustrate the mesh adaptivity, we have included a plot of the mesh trajectories in Figure 6. Each line represents the trajectory of one mesh point in time, and we can see that the mesh points cluster nicely around the edges of the waves as they move.

To illustrate the performance of our methods, we have in Figure 5 compared solutions obtained by using the -preserving moving mesh method with the solutions obtained by using a fourth order Runge–Kutta method on a static mesh, with the same, and quite few, degrees of freedom. The DG2MM solution is visibly closer to the solutions in Figure 4. The non-preserving RK scheme does a worse job of preserving the amplitude and speed of the waves compared to the DG2MM scheme, and we observe unwanted oscillations.

In Figure 7 we have plotted the Hamiltonian errors for this problem. Again we see that the energy preserving schemes preserve both Hamiltonians better than the Runge–Kutta scheme, but we do also observe that the DG1 scheme preserves better than the DG1MM scheme, and vice versa for the DG2 and DG2MM schemes. Note also that an increase in the errors can be observed when the two waves interact, but that this increase is temporary.

## 6 Conclusions

In this paper, we have presented energy preserving schemes for a class of PDEs, first on general fixed meshes, and then on adaptive meshes. These schemes are then applied to the BBM equation, for which discrete schemes preserving two of the Hamiltonians of the problem are explicitly given.

Numerical experiments are performed, using the energy preserving moving mesh schemes on two different BBM problems: a soliton solution, and two waves interacting. Plots of the phase and shape errors illustrate how, for the given parameters, the usage of moving meshes gives improved accuracy, while the integral preservation gives comparable results to existing methods, without yielding a categorical improvement. We will remark, however, that in many cases, the preservation of a quantity such as one of the Hamiltonians in itself may be a desired property of a numerical scheme. For the two wave interaction problem, we do not have an analytical solution to compare to, but plots of the solution indicate that our schemes perform well compared to a Runge–Kutta scheme.

Although the numerical examples presented here are simple one-dimensional problems, the adaptive discrete gradient methods should also be applicable for multi-dimensional problems. This could be an interesting direction for further work, since the advantages of adaptive meshes are typically more evident when increasing the number of dimensions.

## References

- [1] R. Courant, K. Friedrichs, and H. Lewy, “Über die partiellen Differenzengleichungen der mathematischen Physik,” Math. Ann., vol. 100, no. 1, pp. 32–74, 1928.
- [2] S. Li and L. Vu-Quoc, “Finite difference calculus invariant structure of a class of algorithms for the nonlinear Klein-Gordon equation,” SIAM J. Numer. Anal., vol. 32, no. 6, pp. 1839–1875, 1995.
- [3] E. Hairer, C. Lubich, and G. Wanner, Geometric numerical integration, vol. 31 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second ed., 2006. Structure-preserving algorithms for ordinary differential equations.
- [4] O. Gonzalez, “Time integration and discrete Hamiltonian systems,” J. Nonlinear Sci., vol. 6, no. 5, pp. 449–467, 1996.
- [5] D. Furihata and T. Matsuo, Discrete variational derivative method. Chapman & Hall/CRC Numerical Analysis and Scientific Computing, CRC Press, Boca Raton, FL, 2011. A structure-preserving numerical method for partial differential equations.
- [6] E. Celledoni, V. Grimm, R. I. McLachlan, D. I. McLaren, D. O’Neale, B. Owren, and G. R. W. Quispel, “Preserving energy resp. dissipation in numerical PDEs using the “average vector field” method,” J. Comput. Phys., vol. 231, no. 20, pp. 6770–6789, 2012.
- [7] T. Yaguchi, T. Matsuo, and M. Sugihara, “An extension of the discrete variational method to nonuniform grids,” J. Comput. Phys., vol. 229, no. 11, pp. 4382–4423, 2010.
- [8] T. Yaguchi, T. Matsuo, and M. Sugihara, “The discrete variational derivative method based on discrete differential forms,” J. Comput. Phys., vol. 231, no. 10, pp. 3963–3986, 2012.
- [9] Y. Miyatake and T. Matsuo, “A note on the adaptive conservative/dissipative discretization for evolutionary partial differential equations,” J. Comput. Appl. Math., vol. 274, pp. 79–87, 2015.
- [10] S. Eidnes, B. Owren, and T. Ringholm, “Adaptive energy preserving methods for partial differential equations,” Advances in Computational Mathematics, Sep 2017.
- [11] D. Cohen and X. Raynaud, “Geometric finite difference schemes for the generalized hyperelastic-rod wave equation,” J. Comput. Appl. Math., vol. 235, no. 8, pp. 1925–1940, 2011.
- [12] C. Lu, W. Huang, and J. Qiu, “An adaptive moving mesh finite element solution of the regularized long wave equation,” Journal of Scientific Computing, pp. 1–23, 2016.
- [13] R. I. McLachlan, G. R. W. Quispel, and N. Robidoux, “Geometric integration using discrete gradients,” R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci., vol. 357, no. 1754, pp. 1021–1045, 1999.
- [14] W. Huang and R. Russell, Adaptive Moving Mesh Methods, vol. 174 of Springer Series in Applied Mathematical Sciences. Springer-Verlag, New York, 2010.
- [15] I. Babus̆ka and B. Guo, “The h, p and h-p version of the finite element method; basis theory and applications,” Advances in Engineering Software, vol. 15, pp. 159–174, 1992.
- [16] D. Peregrine, “Calculations of the development of an undular bore,” Journal of Fluid Mechanics, vol. 25, no. 02, pp. 321–330, 1966.
- [17] T. B. Benjamin, J. L. Bona, and J. J. Mahony, “Model equations for long waves in nonlinear dispersive systems,” Philos. Trans. Roy. Soc. London Ser. A, vol. 272, no. 1220, pp. 47–78, 1972.
- [18] T.-C. Wang and L.-M. Zhang, “New conservative schemes for regularized long wave equation,” Numerical Mathematics A Journal of Chinese Universities English Series, vol. 15, no. 4, pp. 348–356, 2006.
- [19] J. A. Cottrell, T. J. Hughes, and Y. Bazilevs, Isogeometric Analysis. Wiley, 2009.
- [20] C. J. Budd, W. Huang, and R. D. Russell, “Adaptivity with moving grids,” Acta Numer., vol. 18, pp. 111–241, 2009.
- [21] J. Blom and J. Verwer, “On the use of the arclength and curvature monitor in a moving-grid method which is based on the method of lines,” tech. rep., NM-N8902, CWI, Amsterdam, 1989.
- [22] W. Craig, P. Guyenne, J. Hammack, D. Henderson, and C. Sulem, “Solitary water wave interactions,” Physics of Fluids, vol. 18, no. 5, p. 057106, 2006.