On time relaxed schemes and formulations for dispersive wave equations

1.   Introduction

Nonlinear dispersive waves arise in various fields of science including the solid ^[29] and fluid mechanics ^[20]. Below we focus on models stemming from the free surface hydrodynamics mainly due to scientific interests of the authors. In this field, analytical solutions (such as solitary or cnoidal waves ^[8,10]) are available only for some simplified models. That is why the numerical simulation remains one of the main tools in nonlinear dispersive wave studies. In this work the authors have an ambition to propose a novel strategy to tackle numerically these problems.

Let us review first the modern numerical approaches to this problem. Historically, the finite differences were applied to dispersive wave problems ^[19]. Today, the trend is to apply GALERKIN–type methods (e.g. FEM) in smooth situations ^[1,26] and finite volumes ^[14,15] or coupled finite volumes (for the hyperbolic part)/finite differences (for the dispersive part) via the operator splitting ^[6]. Notice that there is an effort to develop discontinuous GALERKIN–type methods for dispersive wave equations as well ^[16,36], which combine the robustness of finite volumes with the accuracy of finite elements. However, the 'price to pay' (i.e. the CPU time) remains excessively high turning these methods into a luxury limousine rather than a working horse of numerical analysis. Finally, for periodic situations one can apply highly efficient FOURIER-type pseudo-spectral methods ^[13]. On the borderline between finite differences and pseudo-spectral methods there exist compact finite difference schemes proposed by LELE (1992) with spectral-like resolution. To give an example, these schemes were applied to the SERRE equations in ^[9].

In the present study we propose to change the numerical strategy in contrast to studies described above. The main idea consists in modifying the governing equations by slightly perturbing them with some ad-hoc terms. In this way, we gain the structure suitable for numerical simulations and we hope that solutions of the perturbed and un-perturbed problems will remain close provided that perturbation parameter is small. This idea is inspired by pseudo-compressibility methods proposed for the numerical simulation of incompressible NAVIER–STOKES equations (see e.g. ^[21]). The main philosophy of this approach is as follows. When we solve numerically a differential equation, some errors are introduced due to the underlying discretization process. So, perhaps, for the sake of convenience, one can perturb also the governing equations within the same discretization error which is already introduced. And if everything is done judiciously, the end user will not even see the difference between the relaxed and the original problems solutions. It makes a lot of 'ifs', but nevertheless we undertake this programme below for some well-known dispersive wave equations. We would like to mention here also another relaxation technique based on the application of local time-averaging operators ^[2]. A relaxation scheme for the KdV equation was proposed in ^[4]. However, their formulation contains 2^nd order derivatives. Below, we will propose an alternative formulation which has an advantage to involve only the 1^st order derivatives. The idea we employ can be traced back at least to various local^* formulations used in continuous ^[33,34] and discontinuous ^[24,36] GALERKIN methods. However, we push it one step further towards the so-called relaxed local formulations.

^* In the case of continuous GALERKIN methods this formulation is also sometimes called the modified one.

The present study is organized as follows. In Section 2 we propose the relaxed formulations for several well-known dispersive wave models. Some mathematical properties of a relaxed formulation for the KdV equation are discussed in Section 3. Our numerical approach based on this relaxed formulation of the KdV equation is outlined in Section 4. The first numerical results are presented in the same Section. Finally, in Section 5 we outline the main conclusions and perspectives of our study. The relaxation of some other dispersive wave equations is discussed in Appendix A.

2.   Mathematical models

Instead of working out the most general situation, in the present study we follow the so-called GELFAND^† principle which states that a theory should be illustrated on the simplest non-trivial example. Below we consider one such example and two others are treated in Appendix A. All cases which steem from the water wave theory ^[32]. For each model we develop the corresponding relaxed formulation.

^† Here we mean I. M. GELFAND to avoid any confusion.

2.1.   Korteweg–de Vries equation

The celebrated KORTEWEG–DE VRIES (KdV) equation was derived independently by J. BOUSSINESQ ^[7] and D. KORTEWEG & G. DE VRIES ^[22] at the end of the $XIX^{\, {\rm{th}}}$ century and in the scaled form it reads:

The variable $u\, (x, \, t)$ may be interpreted physically as the free surface elevation or the horizontal depth-integrated fluid velocity. There exist transformations in order to obtain the canonical form (4.1) ^[20]. We shall rewrite the last equation in the conservative form in order to obtain the flux of the evolution variable $u\,$:

Now, the order of derivatives has to be lowered. It can be done similarly to the ODE theory by introducing additional variables:

At this step we are compatible with local formulations used in the framework of discontinuous GALERKIN methods ^[24,36]. It is obvious that the last system is completely equivalent to equation (2.2) (at least for the smooth solutions). However, the last system is not of the evolution type. It is here that the relaxation comes into the play in order to correct this shortcoming:

where $\delta\ \ll\ 1$ is a small parameter. It can be seen that in the limit $\delta\ \to\ 0$ we recover the original formulation (2.1) (at least formally).

Equations (2.4), (2.5) allow also to assess the physical meaning of the parameter $\delta\,$. Indeed, the terms $v$ and $\delta\, v_{\, t}$ (and $w$ with $\delta\, w_{\, t}$) need to have the same units. Thus, the multiplication by $\delta$ has to compensate the derivative with respect to time. Henceforth, from dimensional arguments, $\delta$ has the unit of time.

3.   Properties

Above we presented relaxed formulations for three well-known equations. However, the properties of new formulations have to be studied. From now on we focus on the KdV equation only (2.3)–(2.5) in order to understand the relaxation effects exhaustively. Other equations will be tackled in following studies.

3.1.   Dispersion relation analysis

In order to study the dispersive properties, first we have to linearize equations (2.3)–(2.5):

Now, we look for plane-wave solutions of the form

where $k$ is the wavenumber, $\omega$ is the frequency and $\bigl\{u_{\, 0}, \, v_{\, 0}, \, w_{\, 0}\bigr\}\ \in\ \mathbb{R}$ are some real amplitudes. By substituting^‡ the plane-wave ansatzä into a linear system (3.1)–(3.3), we obtain a system of linear algebraic equations, where unknowns are $\bigl\{u_{\, 0}, \, v_{\, 0}, \, w_{\, 0}\bigr\}\,$:

^‡ This operation can be seen as the FOURIER transform in both space and time.

The necessary condition to have non-trivial solutions to the above linear system of equations is that its determinant is equal to zero:

After some simplifications we obtain the desired dispersion relation between $\omega$ and $k\,$:

So, one can see that by taking the limit $\delta\ \to\ 0$ in the last equation, we recover the classical KdV dispersion relation $\omega\ = \ -\, k^{\, 3}\,$. This polynomial equation (3.4) in wave frequency $\omega$ admits in general three complex roots. Their general expression is rather cumbersome. However, we can construct their asymptotic expansions:

These expansions (3.5), (3.6) might be used to interpret the behaviour of numerical solutions. Since we deal with singular perturbations in equation (3.4), the passage to the limit $\delta\ \to\ 0$ in solutions $\omega_{\, 2, \, 3}\, (k)$ is far from being trivial. However, from the expansion (3.5) we learn that parameter $\delta$ plays the both rôles of dissipation and dispersion since it enters into real and imaginary parts of the regular branch $\omega_{\, 1}\, (k)\,$. We would like to highlight the fact that the relaxation parameter $\delta$ affects the dispersive properties of the regular branch (i.e. $\operatorname{Re}\omega_{\, 1}\, (k)$) only to the second order in $\delta\,$. Thus, we can conclude that the proposed relaxation approach is very careful with respect to the dispersion relation of the original KdV equation.

Remark 1. Taking the limit in the polynomial equation (3.4) turns out to be an easy task, as it often happens, than taking the limit $\delta\ \to\ 0$ in the solutions $\omega_{\, 2, \, 3}\, (k)\,$. On the other hand, the branch $\omega_{\, 1}\, (k)$ seems to be completely regular from this point of view. Two other branches $\omega_{\, 2, \, 3}\, (k)$ are artificial modes introduced by relaxation.

4.   Numerical methods and results

In this Section we are going to study numerically the relaxed formulations. For this purpose we propose three different discretisations of the KdV equation. One of them is based on the relaxed formulation proposed above and two others are popular finite difference schemes (CRANK–NICOLSON and SANZ–SERNA) based on the original KdV equation. The numerical advantages of the relaxed formulation are illustrated below.

Before considering the time schemes, we focus on the discretization in space which is realized with finite differences. The use of compact schemes as presented in ^[23] allows to reach a high order of accuracy, with a spectral-like resolution, while preserving the capabilities of the finite differences approach, including the implementation on a general CARTESIAN grids and for a variety of boundary conditions. When dealing with surface wave equations, the high level of accuracy is an important property to capture correctly the active frequencies of the solution without working on very fine grids. Hence, the computation can be done with a reasonable computational time.

At first we recall briefly the principle of the Compact Schemes (CS) and we give the discretization matrices that will be used for the simulation of the models presented in Section 2.

In two words, the compact schemes consist in approaching a linear operator (differentiation as well as interpolation) by a rational (instead of polynomial-like) finite differences scheme. These schemes are then implicit, this allows to increase the accuracy and mimic the spectral global dependence, see ^[23] for more details. The approximations at grid points of the operators applied to a regular function $u$ is realized as follows.

Let $U\ \mathop{\stackrel{\, \mathit{\rm{def}}}{: = }\, }\ (U_{\, 1}, \, \ldots, \, U_{\, N})^{\, \top}$ be a vector whose the components are the approximations $u$ at (regularly spaced) grid points $x_{\, i}\ = \ i\, h\,$, $i\ = \ 1, \, \ldots, \, N\,$, here $h\ = \ \dfrac{1}{N}$ is the spatial step-size. We compute approximations of $V_{\, i}\ = \ \mathcal{L}\, (u)\, (x_{\, i})$ as solution of a system

the approximation matrix is then formally $B\ = \ P^{\, -1}\cdot Q\,$. For simplicity we present here fourth order accurate CS. When considering periodic boundary conditions, the matrices $P$ and $Q$ are for the first order derivative in space:

For the second order derivative in space we have the following pair of matrices:

Finally, for the third order derivative in space we have

We now present the semi-implicit relaxed time scheme we will use for the simulation, stressing out its enhanced stability as compared to classical semi-implicit schemes (CRANK–NICOLSON's for the linear terms and forward EULER's for the nonlinear ones) while giving comparable results to the ones obtained by fully nonlinear SANZ–SERNA's which is second order accurate in time and unconditionally stable.

Let $\mathcal{D}_{\, x}\,$, $\mathcal{D}_{\, xx}$ and $\mathcal{D}_{\, xxx}$ be the discretization matrices of the first, second and third order derivative on a $N$ regularly spaced grid points, with periodic boundary conditions; these matrices will be obtained in practice with compact schemes as presented above; $\mathbb{I}_{\, N}$ will denote the $N\times N$ identity matrix. The relaxed scheme produces the iterations

We rewrite this coupled system by blocs as

We can now introduce the iteration matrix

We can resume the relaxed scheme as follows

We now present the two references schemes to which we will compare the relaxed one. We set

and we define the corresponding iteration matrix $\mathcal{M}_{\, \mathit{\rm{CN}}}$ as

We will use also

4.1.   Linear stability analysis

Before comparing the stability properties of the classical CN and of the relaxed schemes, we give hereafter a simple but instructive result illustrating the advantage of the compact schemes in a context in which the spectral properties of the operators must be restituted at the discrete level. Writing the CN scheme for the AIRY equation, we have the following induction relation among FOURIER coefficients attached to the $m$^th frequency

hence

The FOURIER symbol of the CN operator is $\dfrac{1\ +\ \mathit{\rm{i}}\; \dfrac{\Delta t}{2}\; m^{\, 3}}{1\ -\ \mathit{\rm{i}}\; \dfrac{\Delta t}{2}\; m^{\, 3}}$ and its values are all displayed on the unit circle, for any $\Delta t\ > \ 0\,$. It is then an important property to be captured by the spectrum of the CRANK–NICOLSON matrix, when dealing with finite differences. We give hereafter in Figure 1 the spectra of the CN matrix when the 2^nd order FD and the 4^th order compact FD are used. We observe that the compact scheme allows to capture numerically the spectral properties of the linear operator as in the analytical FOURIER analysis, while the 2^nd order finite differences succeed to capture correctly only a small part of the spectrum. Comparable results are obtained in Figure 2 when considering 6^th order compact FD: they illustrate the capabilities of the compact schemes to mimmic the spectral properties of the linear propagation operator.

4.2.   Enhanced stability analysis

A first way to understand the effect of the relaxation in the stability of the time marching schemes is to consider the linear system (3.1). We have used the 6^th order compact schemes for the discretization of the spatial derivatives, as presented in the beginning of Section 4. We compare here in the complex plane the eigenvalues of the two matrices used in the numerical schemes: $\mathcal{M}_{\, r}$ and $\mathcal{M}_{\, \mathit{\rm{CN}}}\,$. We recall that these matrices are not of the same size so there is no reason that their relative spectra coincide. The eigenvalues of $\mathcal{M}_{\, \mathit{\rm{CN}}}$ are mapped on the unit circle (see Figure 1), we present hereafter comparison with those of $\mathcal{M}_{\, r}$ for different values of $N\,$, $\Delta t$ and $\delta$ (the relaxation parameter). At first, we take a very small value of $\delta\ = \ 10^{\, -17}\,$, however the results we obtain are identical all for larger but still small values of $\delta\,$, up to $10^{\, -5}\,$, as illustrated in Figure 3. For larger values of $\delta\,$, say $\delta\ = \ 10^{\, -4}\,$, $\delta\ = \ 10^{\, -3}\,$, $\delta\ = \ 10^{\, -1}\,$, eigenvalues of $\mathcal{M}_{\, r}$ are placed outside the unit disk making the relaxed scheme unstable, see Figures 4, 5 and 6.

Below, in Table 1, we give the maximum and the minimum of the modulus of the eigenvalues of the $\mathcal{M}_{\, r}$ and $\mathcal{M}_{\, \mathit{\rm{CN}}}$ matrices for different values of the time step $\Delta t$ and relaxation parameter $\delta\,$. We can observe the influence of $\delta$ on the spectral radius of $\mathcal{M}_{\, r}\;$: the relaxed scheme becomes unstable, i.e. $\rho\, (\mathcal{M}_{\, r})\ > \ 1\,$, for not sufficiently small values of $\delta\,$, say, e.g. $\delta\ \geqslant\ 10^{\, -4}\,$; for small values of $\delta\,$, we have $\rho\, (\mathcal{M}_{\, r})\ = \ 1\,$, the eigenvalues of $\mathcal{M}_{\, r}$ are perfectly matched on the unit circle and the scheme is then unconditionally stable.

4.3.   Simulation of a soliton

To illustrate the effect of the time relaxation, we consider the simulation of a soliton on the domain $I\ = \ \bigl[\, 0, \, \ell\, \bigr]\,$, for different values of $\delta$ and of $\Delta t\,$, namely

with $c\ = \ \frac{a}{3}$ and $\kappa\ = \ \sqrt{c}\,$. In our experiments presented below we use $a\ = \ 0.8\,$, $\ell\ = \ 100$ and $x_{\, 0}\ = \ \dfrac{\ell}{2}\,$. The interval $I$ is discretized with $N\ = \ 400$ equidistant points and the space discretization is realized by using 6^th compact schemes. Of course, due to the nonlinear orbital stability properties of the solitons ^[35], the approximation to the solitary wave by the numerical solution of the relaxed system cannot be considered in a long time interval, making the validation delicate. Indeed, the classical error measures such as $L_{\, p}$ norms are irrelevant for orbits, meanwhile a numerical method approximates a classical solution and not an orbit.

We observe that the relaxation allows to approach the exact solution with an expected level of accuracy. For instance, when $\delta\ = \ 10^{\, -\, 17}\,$, the solutions computed by the different schemes coincide on a fairly long time interval. With $\delta\ = \ 10^{\, -\, 4}$ and $\Delta t\ = \ 0.1\,$, the solution coincide on a smaller time interval as it is expected. More generally, we notice in our numerical experiments that the pairs $(\delta\, ,\, \Delta t)$ that make the relaxed scheme stable in the linear case are operant for the nonlinear KdV equation as well^§. The contrary is also observed, a numerical instability holds, for example, for $\Delta t\ = \ 0.005$ and $\delta\ = \ 0.001\,$, see Table 1. The numerical result for the soliton propagation are shown in Figure 7 for two values of the relaxation parameter $\delta\ = \ 10^{\, -\, 17}$ and $\delta\ = \ 10^{\, -\, 4}\,$. On the right panels we show the evolution of the $L_{\, 2}$ error computed thanks to the exact Solution (4.1).

^§ However, we do not have a mathematical proof for this statement.

4.4.   Second order relaxed times schemes

As illustrates above, relaxed schemes allow a correct numerical integration using only approximation to $\partial_{\, x}$ which is generally easily available in a number of situations, while the approximation to $\partial_{\, x\, x\, x}$ can be very tricky to build. However, the relaxation scheme is only first order accurate in time and as the semi-implicit EULER scheme, it reveals not to be suited for long time interval simulations. To overcome this drawback we propose here to reach a second order of accuracy by using a classical RICHARDSON extrapolation. In two words, the numerical time integration of ODE

by the forward EULER scheme defines the iterations

which are first order accurate approximations to $u\, (k\, \Delta t)\,$. The RICHARDSON extrapolated sequence is defined by

and is second order accurate in time. We will start here from a simple IMEX method, says backward EULER for the linear terms and forward for the nonlinear ones: in other words if $F\, (u)$ writes as $F\, (u)\ = \ -\, A\, u\ +\ H\, (u)\,$, then the propagator is formally defined as

The IMEX relaxed scheme consists then in solving the following linear system at each step:

The matrix of the system is noted by $M^{\, (\mathit{\rm{IMEX}})}_{\, \Delta t\, ,\, \delta}\,$. The Extrapolated Relaxed schemes writes as

Hereafter, in Figures 8–9, we present the comparison of the evolution of the soliton and its numerical approximations (SANZ–SERNA and extrapolated Relaxed and IMEX schemes). We observe that the time extrapolation allows to approach the exact solution with a good accuracy on longer time intervals, when considering different values of $\delta$ and $\Delta t\,$; this is notable particularly, e.g. when taking $\delta\ = \ 10^{\, -4}$ and $\Delta t\ = \ 10^{\, -1}\,$.

5.   Discussion

Above we presented some rationale behind time-relaxed formulations for several well-known dispersive wave equations. The main conclusions and perspectives of our study are outlined below.

5.1.   Conclusions

In this study we presented an approximate reformulation for several dispersive wave equations. This formulation was inspired somehow by quasi– (or pseudo–) compressibility methods to solve incompressible NAVIER–STOKES equations ^[21]. So, by following this 'philosophy' we proposed relaxed formulations for celebrated KORTEWEG–DE VRIES (KdV), BENJAMIN–BONA–MAHONY (BBM) and PEREGRINE's system of equations. However, it is obvious that the same technique can be extended to many other scalar and vectorial dispersive wave equations. This formulation has a simple advantage to involve first order derivatives only while being in the form of a coupled evolution problem. This is the main difference with various local (or modified) reformulations used in continuous and discontinuous GALERKIN methods ^{[24,33,34,36]}. Many standard numerical methods can be applied to solve the proposed relaxed formulation numerically. In the present study we applied compact finite differences (with spectral-like resolution) to discretize the problem in space along with a simple time stepping. The presented numerical tests and validations for the KdV equation show that the relaxed discrete formulation possesses a larger CFL–type stability limit comparing to similar compact discretizations of the standard KdV equation (without relaxation). This preliminary conclusion indicates that relaxed schemes might be good candidates for the numerical simulation of notoriously stiff systems of equations such as the KdV–KdV BOUSSINESQ-type system extensively studied numerically in e.g. ^[5].

5.2.   Perspectives

Above we presented some numerical illustrations for the classical KdV equation only. As it was mentioned in the previous Section, this formulation was extended to other weakly-nonlinear models as well. The main goal consists in extending these time-relaxation numerical methods to BOUSSINESQ-type systems of weakly nonlinear weakly dispersive equations such as the celebrated classical PEREGRINE ^[28], modified PEREGRINE ^[12] or even fully nonlinear SERRE–GREEN–NAGHDI equations ^[31] (some more conventional numerical strategies for these equations were outlined in e.g. ^[13,25,26]).

On a more theoretical side, we would like to obtain the estimations of the difference between the perturbed and unperturbed problems solutions. In other words, we would like to have theoretical arguments to state that the proposed relaxation process generates solutions which remain close to those of the original equation. This is the main point on our 'theoretical' agenda.

Acknowledgments

D. DUTYKH would like to acknowledge the hospitality of the Laboratory LAMFA and of the University of PICARDIE JULES VERNE during his visit in November 2016. Reciprocally, J.-P. CHEHAB acknowledges the hospitality of the Laboratory of Mathematics (LAMA UMR #5127) and of the University SAVOIE MONT BLANC during his visit in December 2017.

Conflict of interest

The authors declare that there is no conflicts of interest in this paper.

A.   Other types of dispersive equations

In this Appendix we provide two other examples of relaxation for some widely used dispersive wave equations (scalar and system case).

A.1.   Benjamin–Bona–Mahony equation

The celebrated BENJAMIN–BONA–MAHONY (BBM) equation was derived first in ^[3,27] and can be recast in the following dimensionless form:

Similarly to the KdV case the variable $u\, (x, \, t)$ can be the free surface elevation or a horizontal fluid velocity either. It is possible to propose a similar relaxed formulation for the BBM equation as well. Indeed, let us recast it first in a conservative form:

The 2^nd derivative can be lowered by introducing additional variables:

It is convenient to introduce a new variable $\beta \, (x, \, t)\ \mathop{\stackrel{\, \mathit{\rm{def}}}{: = }\, }\ \bigl(u\ -\ w\bigr)\, (x, \, t)$ to simplify the first equation:

The last step consists in adding relaxation terms to obtain an evolutionary system in all variables:

with $\delta\ \ll\ 1$ being again a small parameter.

A.2.   Peregrine system

The PEREGRINE system was proposed by D. H. PEREGRINE (1967) in ^[28]. In a particular case of a fluid layer of constant depth (i.e. the even bottom case) this system reads

where $d\ > \ 0$ is the constant fluid depth and $g\ > \ 0$ is the gravity acceleration. The variables $\eta\, (x, \, t)$ and $u\, (x, \, t)$ are the free surface elevation and the depth-averaged horizontal velocity correspondingly. The sketch of the fluid domain is shown in Figure 10.

As the first step, the PEREGRINE system (A.1), (A.2) is rewritten in the following conservative form:

The last conservative PEREGRINE system can be lowered and relaxed similar to the BBM case (see the Section A.1) by introducing new variables:

where $\delta\ \ll\ 1$ and $\rho\, (x, \, t)\ \mathop{\stackrel{\, \mathit{\rm{def}}}{: = }\, }\ u\, (x, \, t)\ -\ { {1\over3}}\, d^{\, 2} w\, (x, \, t)$ is an auxiliary variable. If one takes the limit $\delta\ \to\ 0\,$, we shall recover a formulation similar to what is used in local (or modified) continuous GALERKIN methods ^[33,34].

Remark 2. The celebrated KdV and BBM models can be written for the variable $u\, (x, \, t)$ being the free surface elevation $\eta\, (x, \, t)$ or the horizontal velocity (after a few changes of variables). In this way these equations appear in the theory of water waves. However, the KdV and BBM-type equations appear in many other physical settings as well (see e.g. ^{[11,17,18,20,30]}).