# Numerical solutions to a BBM-Burgers model with a nonlocal viscous term

Serge Dumont, Imen Manoubi

► To cite this version:

Serge Dumont, Imen Manoubi. Numerical solutions to a BBM-Burgers model with a nonlocal viscous term. Numerical Methods for Partial Differential Equations, 2018, 34 (6), pp.2279-2300. 10.1002/num.22291 . hal-01818121

HAL Id: hal-01818121

https://hal.science/hal-01818121

Submitted on 13 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - NonCommercial 4.0 International License
# Numerical solutions to a BBM-Burgers model with a nonlocal viscous term

Serge Dumont<fnref n="1" />,<fnref n="2" /> | Imen Manoubi<fnref n="3" />

## Correspondence

Imen Manoubi, Unité de recherche: Multifractales et Ondelettes, Faculté des Sciences de Monastir, Av. de l'environnement 5019 Monastir, Tunisia. Email: imen.manoubi@yahoo.fr

## Funding information

This work was partially supported by a grant from the Simons Foundation.

In this paper, we numerically investigate the BBM-Burgers equation with a nonlocal viscous term

$$u_t + u_x - \beta u_{txx} + \frac{\sqrt{v}}{\sqrt{\pi}} \frac{\partial}{\partial t} \int_0^t \frac{u(s)}{\sqrt{t-s}} ds + \gamma u u_x = \alpha u_{xx},$$

where $\frac{1}{\sqrt{\pi}} \frac{\partial}{\partial t} \int_0^t \frac{u(s)}{\sqrt{t-s}} ds$ is the Riemann-Liouville half derivative. In particular, we implement different numerical schemes to approximate the solution and its asymptotical behavior. Also, we compare our numerical results with those given in [1, 2] for similar models.

## KEYWORDS

BBM-Burgers equation, decay rate, fractional derivatives, Gear scheme, nonlocal viscous model, quadrature methods, water waves

# 1 | INTRODUCTION

The mathematical modeling and analysis of water wave propagation are challenging topics. In their work, J. Bona et al. have derived a family of Boussinesq systems from the two-dimensional Euler equations for free-surface flow in [3]. Modeling the effects of viscosity on the propagation of long waves is an important challenge that has been investigated since the time of Stokes and has received a lot of interest in the last decade (see [4, 5] and references therein). Besides, P. Liu and T. Orfila [6], D. Dutykh, and F. Dias [7] have independently derived viscous asymptotic models for transient long-wave propagation including viscous effects. These effects appear as nonlocal terms in the form of convolution integrals. The derivation of this model holds in 3 D and 2 D cases. Using a one-way wave reduction (see [3, 8] for details), the authors in [9] investigated a reduced nonlinear model that reads

<footnote n="1">¹University of Nîmes, Pl. Gabriel Péri, 30000 Nîmes, France</footnote>

<footnote n="2">²IMAG, University of Montpellier, CNRS UMR, 5149 Pl. E. Bataillon, 34095, Montpellier cedex 5, France</footnote>

<footnote n="3">³Unité de recherche: Multifractales et Ondelettes, Faculté des Sciences de Monastir, Av. de l'environnement 5019 Monastir, Tunisia</footnote>
$$u_t + u_x + \beta u_{xxx} + \frac{\sqrt{v}}{\sqrt{\pi}} \int_0^t \frac{u_t(s)}{\sqrt{t-s}} ds + \gamma u u_x = \alpha u_{xx}, \quad (1)$$

where $\frac{1}{\sqrt{\pi}} \int_0^t \frac{u_t(s)}{\sqrt{t-s}} ds$ is the Caputo half-derivative. Here $u$ is the horizontal velocity of the fluid, $-\alpha u_{xx}$ is the usual diffusion, $\beta u_{xxx}$ is the geometric dispersion, $\frac{\sqrt{v}}{\sqrt{\pi}} \int_0^t \frac{u_t(s)}{\sqrt{t-s}} ds$ stands for the nonlocal diffusive-dispersive term. The parameters $\beta, \nu, \gamma$, and $\alpha$ are dedicated to balance the effects of viscosity and dispersion against nonlinear effects. Moreover, in the recent work [2], one of the authors has considered the following water wave model

$$u_t + u_x + \beta u_{xxx} + \frac{\sqrt{v}}{\sqrt{\pi}} \frac{\partial}{\partial t} \int_0^t \frac{u(s)}{\sqrt{t-s}} ds + \gamma u u_x = \alpha u_{xx}, \quad (2)$$

where $\frac{1}{\sqrt{\pi}} \frac{\partial}{\partial t} \int_0^t \frac{u(s)}{\sqrt{t-s}} ds$ is the Riemann-Liouville half derivative.

Particularly, it is proved the local and the global existence result and decay estimates for the integro-differential equation (2) when $\beta = 0, \nu = \alpha = \gamma = 1$ supplemented with the initial condition $u_0 \in L^1(\mathbb{R}) \cap L^2(\mathbb{R})$. Precisely, the following theorem is stated

**Theorem 1.1** (I. Manoubi, [2]) Let $u_0 \in L^2(\mathbb{R})$, then there exists a unique local solution $u \in C([0, T); L_x^2(\mathbb{R}))$ of (2).

Moreover for $u_0 \in L^1(\mathbb{R}) \cap L^2(\mathbb{R})$, there exists a positive constant $C_0 > 0$ that depends on $u_0$ such that if $\|u_0\|_{L^1(\mathbb{R})}$ is small enough, there exists a unique global solution $u \in C(\mathbb{R}_+; L_x^2(\mathbb{R})) \cap C^{1/2}(\mathbb{R}_+; H_x^{-2}(\mathbb{R}))$ of (2) given by

$$u(t,x) = [K_{RL}(t, \cdot) \star u_0](x) - N \circledast u^2(t,x), \quad (3)$$

where $K_{RL}$ and $N$ are given by

$$K_{RL}(t,x) = \frac{1}{2\sqrt{\pi t}} e^{-\frac{x^2}{4t}} e^{-x} \left( 1 - \frac{1}{2} \int_0^{+\infty} e^{-\frac{\mu^2}{4t} - \frac{\mu|x|}{2t} - \frac{\mu}{2}} d\mu \right),$$

and

$$N(t,x) = \frac{1}{4\sqrt{\pi t}} \partial_x \left( e^{-\frac{x^2}{4t}} e^{-x} \left( 1 - \frac{1}{2} \int_0^{+\infty} e^{-\frac{\mu^2}{4t} - \frac{\mu|x|}{2} - \frac{\mu}{2}} d\mu \right) \right).$$

with $x^- = \frac{|x|-x}{2} = \max(-x, 0)$, $\star$ represents the usual convolution product and $\circledast$ is the time-space convolution product defined by

$$v \circledast w(t,x) = \int_0^t \int_{\mathbb{R}} v(t-s, x-y) w(s, y) ds dy.$$

whenever the integrals make sense. In addition, we have the following estimate

$$\max(t^{1/4}, t^{3/4}) \|u(t, \cdot)\|_{L_x^2(\mathbb{R})} + \max(t^{1/2}, t) \|u(t, \cdot)\|_{L_x^\infty(\mathbb{R})} \le C_0. \quad (4)$$

The proof of this theorem is presented in [2].

In addition, in the recent work [10], the authors succeeded to remove the smallness condition on the initial data in Theorem 1.1. Moreover, they proved the weak convergence to zero of the solution.
Furthermore, in their recent work [1], S. Dumont and J.-B Duval investigated numerically the decay rate for solutions to the following water wave model

$$u_t + u_x - \beta u_{txx} + \frac{\sqrt{\nu}}{\sqrt{\pi}} \int_0^t \frac{u_t(s)}{\sqrt{t-s}} ds + \gamma uu_x = \alpha u_{xx}, \quad (5)$$

The approximation of time fractional operators has received a lot of interest during last decades for their wide application in fluid, in solid mechanics and in visco-elasticity. The formulation of a numerical stable scheme is crucial but also a difficult issue because of the nonlocal feature of such operators. The classical methods used in the literature [11–17] consist in the approximation of these fractional operators using either convolution integrals or the so-called Gear scheme for fractional operators. Recently, number of authors in the automatic community developed an alternative method, called the diffusive realization, which is devoted to *causal* pseudodifferential operators [18–21]. Different applications of this approach can be found in [22–25]. The main idea of this method is to replace the nonlocal operator by a linear differential equation. The resulting diffusive model is infinite dimensional, but local in time. Hence, the new model is more easy to solve for both analytical and numerical points of view.

In this paper, we are interested in the following equivalent Benjamin-Bona-Mahony (BBM) model of (2)

$$u_t + u_x - \beta u_{txx} + \frac{\sqrt{\nu}}{\sqrt{\pi}} \frac{\partial}{\partial t} \int_0^t \frac{u(s)}{\sqrt{t-s}} ds + \gamma uu_x = \alpha u_{xx}. \quad (6)$$

We implement two numerical schemes to approximate the solution of (6). The first one is detailed in [1, 2, 13] and is based on the Gear scheme for the approximation of the Riemann-Liouville half-derivative. The second method is based on the diffusive realization of the nonlocal operator supplemented with a splitting scheme (see [8, 26] and references therein). We perform numerical simulations on the solutions and on the decay rates for different values of the parameters $\beta$, $\nu$, $\gamma$, and $\alpha$. We compare between these schemes. Furthermore, we compare our numerical results with those given in [1, 2, 10].

**Remark 1** We note that the well-posedness of the model (6) may be proved mathematically for initial data $u_0 \in L^2(\mathbb{R})$ using the diffusive realization of the half-order derivative and following the same steps as presented at [10].

The outline of this article is as follows: in Section 2, we develop the dispersion relation of the model (6). Then, in Section 3, we present a first numerical scheme of (6) and numerical results using the Gear scheme to approximate the nonlocal term. In Section 4, we perform a second numerical scheme based on the diffusive realization of the nonlocal term followed by several numerical simulations for model (6). A comparison between the two schemes is also performed.

## 2 | DISPERSION RELATION

We discuss, in this section, the dispersion relation for the linearized asymptotic model. Similarly to [9], we take $\beta = 1$, $\gamma = 0$, $\alpha = \nu$ and we consider a Laplace-Fourier analysis due to the presence of the nonlocal term.

Consider the linear BBM-Burgers equation

$$u_t + u_x - u_{txx} = \nu u_{xx}. \quad (7)$$
We seek a plane wave solution of the form $u(t, x) = v(t)e^{ikx}$ with $v(0) = 0$. Substituting this solution into (7), we get

$$ (1 + k^2)v_t + (\nu k^2 + ik)v = 0. \quad (8) $$

We now apply the Laplace transform to (8). We recall that the Laplace transform in time of a function $v$ of exponential order $\alpha$ is defined by

$$ \mathcal{L}(v)(\tau) = \tilde{v}(s) = \int_{0}^{+\infty} v(t)e^{-\tau t}dt, $$

for all $\tau$ such that $\Re(e(\tau)) > \alpha$. Hence, we get

$$ (1 + k^2)\tau + (\nu k^2 + ik) = 0. \quad (9) $$

The real part of $\tau$ namely, $\Re(\tau) = -\frac{\nu k^2}{1+k^2}$ represents the dissipation relation. The imaginary part, denoted by $\omega = -\operatorname{Im}(\tau) = \frac{k}{1+k^2}$, represents the dispersion relation. In the sequel, we linearize (2). Thus, we get

$$ u_t + u_x - u_{txx} + \frac{\sqrt{\nu}}{\sqrt{\pi}} \frac{\partial}{\partial t} \int_0^t \frac{u(s)}{\sqrt{t-s}} ds = \nu u_{xx}. \quad (10) $$

Substituting the plane wave solution into (10), we get

$$ (1 + k^2)v_t + \frac{\sqrt{\nu}}{\sqrt{\pi}} \frac{\partial}{\partial t} \int_0^t \frac{v(s)}{\sqrt{t-s}} ds + (\nu k^2 + ik)v = 0. \quad (11) $$

We now apply the Laplace transform to (11). Since $v(0) = 0$, we get

$$ (1 + k^2)\tau + \sqrt{\nu\tau} + \nu k^2 + ik = 0. \quad (12) $$

In order to solve equation (12), we consider the change of variables: $\tau = z^2$ using the principal determination of the logarithm, such that $\Re(z) > 0$. Hence

$$ (1 + k^2)z^2 + \sqrt{\nu}z + \nu k^2 + ik = 0. \quad (13) $$

Equation 13 has two solutions. We consider the solution $z$ such that $\Re(z) > 0$, namely

$$ z = \frac{-\sqrt{\nu} + \sqrt{\nu - 4(1 + k^2)(\nu k^2 + ik)}}{2(1 + k^2)} $$

Then

$$ -z^2 = -\tau = \frac{-\nu}{2(1 + k^2)^2} + \frac{\sqrt{\nu}}{2(1 + k^2)^2} \sqrt{\nu - 4(1 + k^2)(\sqrt{\nu} - ik)} + \frac{\nu k^2 + ik}{1 + k^2} $$
By restricting to the regime $v \ll k \ll 1$, we obtain

$$ 
\begin{aligned}
\sqrt{v} - 4(1 + k^2)(vk^2 + ik) &= \sqrt{-4ik}\sqrt{1 + o(1)} \\
&= 2e^{-i \operatorname{sgn}(k)\pi/4} \sqrt{|k|} + o(\sqrt{|k|}).
\end{aligned} 
\quad (14) 
$$

Therefore,

$$ -\tau = \frac{-v}{2(1 + k^2)^2} + \frac{\sqrt{v}}{(1 + k^2)^2} \left( \frac{1}{\sqrt{2}} - i \frac{\operatorname{sgn}(k)}{\sqrt{2}} \right) \sqrt{|k|} \\ 
+ \frac{vk^2 + ik}{1 + k^2} + o \left( \frac{\sqrt{v|k|}}{(1 + k^2)^2} \right). 
\quad (15) $$

This implies

$$ -\tau = \left( \frac{-v}{2(1 + k^2)^2} + \frac{\sqrt{|v|k|}}{\sqrt{2}(1 + k^2)^2} + \frac{vk^2}{1 + k^2} \right) \\ 
+ i \left( \frac{k}{1 + k^2} - \operatorname{sgn}(k) \frac{\sqrt{|v|k|}}{\sqrt{2}(1 + k^2)^2} \right) + o \left( \frac{\sqrt{|v|k|}}{(1 + k^2)^2} \right). 
\quad (16) $$

After simplifications

$$ \operatorname{Im}(-\tau) = \left( k - \operatorname{sgn}(k) \frac{\sqrt{|v|k|}}{\sqrt{2}} \right) + o \left( \frac{\sqrt{|v|k|}}{(1 + k^2)^2} \right). 
\quad (17) $$

We observe that $-\operatorname{Im}(\tau)$ has only one nonlinear term: $-\operatorname{sgn}(k) \frac{\sqrt{|v|k|}}{\sqrt{2}}$ which represents the nonlocal dispersion. This term is coming from the nonlocal viscous effect. We would like to point out, as in [9] for the KdV equivalent model (2), that if $v \ll k \ll 1$ the viscous dispersion is dominant with respect to the geometric dispersion coming from the term $u_{txx}$. Actually, the viscosity provides dissipation that is of importance. Nevertheless, it depends on parameters. We will see in examples below (Table 2) that the term $-u_{txx}$ also plays a role.

## 3 | A FIRST NUMERICAL SCHEME

This section deals with the numerical solution of the nonlinear equation

$$ u_t - \beta u_{txx} + \frac{\sqrt{v}}{\sqrt{\pi}} \frac{\partial}{\partial t} \int_0^t \frac{u(s)}{\sqrt{t-s}} ds = \alpha u_{xx} - u_x - \gamma u u_x, 
\quad (18) $$

supplemented with an initial condition $u_0$. Here, $\beta$, $v$, $\alpha$, and $\gamma$ are non-negative parameters.

### 3.1 | Presentation of the scheme

We develop here a first numerical scheme using the Gear scheme. For that purpose, we follow the approach proposed in [1, 2]. First, we present the outline of the Gear scheme developed by A.-C Galucio et al. in [13]. Let $u$ be a time dependent function known only in its discretized values $u^n$ at each time $t^n$, where $n$ is a positive integer. The function $u^n$ is approximated by $u(t^n)$ with $t^n = n\Delta t$, where $\Delta t$,
<table><caption>**TABLE 1** The first five coefficients $g_{n+1}$ of the formal power series (22)</caption><thead><tr><th>j</th><th>α = 1/3</th><th>α = 1/2</th><th>α = 3/4</th></tr></thead><tbody><tr><td>0</td><td>1</td><td>1</td><td>1</td></tr><tr><td>1</td><td>-<sup>4</sup>&frasl;<sub>9</sub></td><td>-<sup>2</sup>&frasl;<sub>3</sub></td><td>-1</td></tr><tr><td>2</td><td>-<sup>7</sup>&frasl;<sub>81</sub></td><td>-<sup>1</sup>&frasl;<sub>18</sub></td><td><sup>1</sup>&frasl;<sub>12</sub></td></tr><tr><td>3</td><td>-<sup>104</sup>&frasl;<sub>2187</sub></td><td>-<sup>1</sup>&frasl;<sub>27</sub></td><td>-<sup>1</sup>&frasl;<sub>108</sub></td></tr><tr><td>4</td><td>-<sup>643</sup>&frasl;<sub>19683</sub></td><td>-<sup>17</sup>&frasl;<sub>648</sub></td><td>-<sup>1</sup>&frasl;<sub>96</sub></td></tr><tr><td>5</td><td>-<sup>4348</sup>&frasl;<sub>177147</sub></td><td>-<sup>19</sup>&frasl;<sub>972</sub></td><td>-<sup>7</sup>&frasl;<sub>846</sub></td></tr></tbody></table>

which is supposed to be fixed, is the time step. Furthermore, let us introduce a delay operator given by $(\delta^{-}u)^{n} = u^{n-1}$. Let $G$ be the Gear operator defined by

$$G = \frac{1}{\Delta t} \left[ \frac{3}{2} I - 2\delta^{-} + \frac{1}{2} (\delta^{-})^2 \right], \quad (19)$$

that approximates the first derivative of $u$ with respect to time. Then the fractional differential operator $G^{\alpha}$ is given by

$$G^{\alpha} = \frac{1}{\Delta t^{\alpha}} \left(\frac{3}{2}\right)^{\alpha} \left[ I - \frac{4}{3}\delta^{-} + \frac{1}{3}(\delta^{-})^2 \right]^{\alpha},$$

which is directly obtained by evaluating the $\alpha$-power of (19). Then, using Newton binomial formula to compute the term in brackets, we get

$$G^{\alpha} = \frac{1}{\Delta t^{\alpha}} \left(\frac{3}{2}\right)^{\alpha} \sum_{j=0}^{\infty} \sum_{l=0}^{j} \left(\frac{4}{3}\right)^{j} \left(\frac{1}{4}\right)^{l} (-1)^{j} C_{\alpha}^{j} (-1)^{l} C_{j}^{l} (\delta^{-})^{j+l}, \quad (20)$$

where $(-1)^{j}C_{\alpha}^{j}$ is given in terms of the Gamma function

$$(-1)^{j} C_{\alpha}^{j} = \frac{\Gamma(j - \alpha)}{\Gamma(-\alpha)\Gamma(j + 1)}.$$

Then the $\alpha$-derivative of $u$ at each time $t^n$ can be approximated by

$$(G^{\alpha}u)^n = \frac{1}{\Delta t^{\alpha}} \left(\frac{3}{2}\right)^{\alpha} \sum_{j=0}^{\infty} g_{j+1} u^{n-j}, \quad (21)$$

where $g_{j+1}$ are rational numbers. For illustrative purposes, we present the first five $G^{\alpha}$-coefficients in Table 1 for three values of $\alpha$: $\frac{1}{3}$, $\frac{1}{2}$, and $\frac{3}{4}$.

As it was mentioned in [12], the Gear operator (19) leads to a three-level step algorithm, backward in time and second order accurate, for the approximation of classical time derivatives. As a consequence, it is a first-order accurate for the approximation of the half derivative. It is worth to note that numerical tests on the convergence of the Gear scheme have been performed in [13], that confirm this property.
Let us notice that since we consider functions $u$ that are vanishing for $t < 0$, the infinite sum in (21) becomes finite

$$ (G^\alpha u)^n = \frac{1}{\Delta t^\alpha} \left(\frac{3}{2}\right)^\alpha \sum_{j=0}^n g_{j+1} u^{n-j}. \quad (22) $$

To write the numerical scheme associated to (18) we follow [1, 2]. For the first iteration, namely for $n = 0$, we use a Crank-Nicolson discretization in time for the linear terms and a fixed-point method for the nonlinear term. To this end, we rewrite the nonlinear term as $uu_x = \frac{1}{4}(u^2)_x + \frac{1}{4}(u^2)_x$. Then we approximate the first term of the right-hand side explicitly and the second one implicitly. Hence the approximate solution $u^1$ verifies the semidiscret scheme

$$ \begin{aligned} & (1 - \beta \partial_{xx}) \frac{u^1 - u^0}{\Delta t} + \frac{1}{2} \sqrt{\frac{3\nu}{2\Delta t}} ((g_2 + g_1)u^0 + g_1 u^1) \\ & = \alpha \frac{u_{xx}^1 + u_{xx}^0}{2} - \frac{u_x^1 + u_x^0}{2} - \frac{\gamma}{2} \frac{(u^0)_x^2 + (u^1)_x^2}{4}. \end{aligned} \quad (23) $$

Then, for $n \ge 1$, we discretize the right-hand side of (18) as in [27]. We use a Crank-Nicolson discretization in time for linear terms and Adams-Bashforth discretization for the nonlinear term. Hence, the proposed discretized equation in time of (18) reads: for all $n \ge 1$

$$ (1 - \beta \partial_{xx}) \frac{u^{n+1} - u^n}{\Delta t} + \sqrt{\nu} (G^{1/2} u)^n = \alpha \frac{u_{xx}^{n+1} + u_{xx}^n}{2} - \frac{u_x^{n+1} + u_x^n}{2} - \frac{\gamma}{2} \frac{(3u^n)_x^2 - (u^{n-1})_x^2}{4}, \quad (24) $$

where

$$ \begin{aligned} (G^{1/2} u)^n &= \frac{1}{2} G^{1/2} (u^{n+1} + u^n) \\ &= \frac{1}{2} \sqrt{\frac{3}{2\Delta t}} \left( \sum_{j=0}^{n+1} g_{n+2-j} u^j + \sum_{j=0}^{n} g_{n+1-j} u^j \right). \end{aligned} $$

In the case $\nu = 0$, this scheme has local truncation error of order $(\Delta t)^2$ and a second-order convergence is observed (for more details, see [27]).

Applying the Fourier transform in space to (23)–(24) provides the scheme:

$$ \begin{aligned} & (1 + \beta \xi^2)(\hat{u}^1 - \hat{u}^0) + \frac{1}{2} \sqrt{\frac{3\nu\Delta t}{2}} ((g_2 + g_1)\hat{u}^0 + g_1 \hat{u}^1) \\ & = -\frac{\Delta t(\alpha\xi^2 + i\xi)}{2}(\hat{u}^1 + \hat{u}^0) - \frac{i\gamma\Delta t\xi}{8}((\hat{u}^0)^2 + (\hat{u}^1)^2). \end{aligned} \quad (25) $$

for all $n \ge 1$

$$ \begin{aligned} & (1 + \beta \xi^2)(\hat{u}^{n+1} - \hat{u}^n) + \frac{1}{2} \sqrt{\frac{3\nu\Delta t}{2}} \left( \sum_{j=0}^{n+1} g_{n+2-j} \hat{u}^j + \sum_{j=0}^{n} g_{n+1-j} \hat{u}^j \right) \\ & = -\frac{\Delta t(\alpha\xi^2 + i\xi)}{2}(\hat{u}^{n+1} + \hat{u}^n) - \frac{i\gamma\Delta t\xi}{8}(3(\hat{u}^n)^2 - (\hat{u}^{n-1})^2). \end{aligned} \quad (26) $$

We note that this first numerical scheme is of order 1 in time. Moreover, this scheme has local truncation error of order $\Delta t$ (see [7, 28], so that first order convergence is expected.
<figure><img src="image_2.png" /><figcaption>FIGURE 1 Error of the time discretization using the Gear scheme [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

## 3.2 | Numerical results and discussion

Similarly to [1], we choose an initial datum $u_0$ which provides an exact BBM soliton for $\alpha = \nu = 0$ and $\beta = \gamma = 1$.

$$u_0(x) = 3(p-1)\mathrm{sech}^2\left(\frac{1}{2}\sqrt{\frac{p-1}{p}}(x-x_0)\right). \quad (27)$$

We take $x_0 = 100$ and $p = 2$. For the numerics, we consider periodic boundary conditions in space over an interval $[0, L]$ with $L$ large enough. In all the simulations hereinafter, we take $L = 800$, the space step size $h = 0.1$ and the time step size $\Delta t = 0.1$. We note that we expect the decay rate of the solution to be as $\|u(t, \cdot)\|_{L^2} \approx Ct^a$ or $\|u(t, \cdot)\|_{L^\infty} \approx Ct^a$ for $t$ large with $a, a' < 0$. Thus, we have the following estimates on the ratios

$$\lim_{t \to \infty} R_2 = \lim_{t \to \infty} \log \left( \frac{\|u(t + \Delta t, \cdot)\|_{L^2}}{\|u(t, \cdot)\|_{L^2}} \right) \left( \log \left( \frac{t + \Delta t}{t} \right) \right)^{-1} = a,$$

$$\lim_{t \to \infty} R_{\infty} = \lim_{t \to \infty} \log \left( \frac{\|u(t + \Delta t, \cdot)\|_{L^{\infty}}}{\|u(t, \cdot)\|_{L^{\infty}}} \right) \left( \log \left( \frac{t + \Delta t}{t} \right) \right)^{-1} = a'$$

In Figure 1, we justify the convergence in time of the numerical scheme (25)–(26). To this end, we fix the parameters values to $\alpha = \beta = 1$, $\gamma = 0.5$, and $\nu = 0.1$. Since we do not know the analytical solution, we denote by $u_{\mathrm{Ref}}^n$ the reference solution when $\Delta t = 0.04$. Let $u^n$ be the numerical solution when increasing the time step $\Delta t$ from 0.05 to 0.2. Then, we denote by $E^n(\Delta t) = \|u_{\mathrm{Ref}}^n - u^n\|_2$ the $L^2$-norm of the error in time. We recall that the solutions are calculated up to time $T = 100$. We plot $E^n$ with respect to $\Delta t$. We observe that the error $E^n$ is decreasing when the time step $\Delta t$ is decreasing. Besides, this figure deals with the order in time of the scheme (25)–(26) which is given by the slope of the curve. We see that the measured values are close to a straight line with slope 1.5. This means that the convergence with respect to the time step of discretization is faster for this example than the expected one, equal to $O(\Delta t)$.
In Figure 2, we simulate the scheme (25)–(26) with parameters values $v = 1$, $\beta, \alpha$ and $\gamma = 0$ or 1. We present the solution at time $T = 500$ and the ratios $R_2$ and $R_\infty$. We observe that the shapes of the numerical solutions are very close. We conclude that the influence of the parameters $\alpha, \gamma$, and $\beta$ is less significant than the nonlocal term on the values of the solution. However, these parameters influence the decay rates. Moreover, in the case $\beta = 0$, the decay rates of the numerical solution match the theoretical results, presented in [2], very well. In addition, when comparing with the equivalent Caputo model (5) investigated in [1], we see that in our case the solution decreases significantly. In fact, in Figure 7 of [1], the maximum values of the numerical solutions are between 0.35 and 0.45. However, in our case, the maximum values of the numerical solutions are between 0.01 and 0.015. A similar observation was done when considering the K.d.V-like equation with the Caputo and the Riemann-Liouville nonlocal terms. For more details, see [2, 9]. Also, we observe that the velocity of the wave (defined e.g., by the evolution of the maximum of the solution) when using Caputo term is greater than that with the Riemann-Liouville term. In fact, in Figure 7 of [1], we observe that the solutions with Caputo term are centered at $x = 255$. However, in our case, the solutions are centered at $x = 130$.

In Figure 3, we simulate the solutions and the ratios $R_2$ and $R_\infty$ when $v = 0$ for different values of the parameters $\beta, \alpha$, and $\gamma$. We observe that the wave moves faster than the case $v = 1$ (see Figure 4). Moreover, the amplitude of the solutions in Figure 3 is greater than in the case $v = 1$ in Figure 4. Also, we observe that in the absence of the nonlocal term, the parameter $\gamma$ describing the nonlinear term plays an important role in this simulation.

However, the parameter $\beta$ affects weakly the amplitude of the solution. In addition, numerical results of decay rates in Figure 3 match very well the theoretical ones for this case. For more details, we refer the reader to [9, 29] and references therein.

In Figure 4, we study the influence of different parameters on the solution and on the decay rates in $L^2$ and in $L^\infty$ norms. We observe that when the $v = 0.1$, the wave moves faster than the case $v = 1$. Similarly to the Figure 3, the parameter $\beta$ does not play an important role in this simulation. In addition, comparing to the results given in [1], we observe that the solution with the Caputo half-derivative moves faster than that with the Riemann-Liouville half-derivative.

In Table 2, we display the values of the decay rates in $L^2$ and in $L^\infty$ norms for different values of the parameters when $\Delta t = h = 0.1$. We observe that when $\beta = 0$, the numerical results match well the mathematical results established in [2] for Equation 2. In addition, when $\beta \neq 0$ the decay rate in $L^2$ norm is about -0.75 and is around -1 for the $L^\infty$-norm. Hence, we deduce that the decay rates of (18) is close to that of (2) but it is different. This difference is due to the dispersion term.

Finally, we calculate the computational time elapsed to simulate the numerical scheme when varying the parameters values. Results are displayed in Table 2. We see that using the Gear scheme is relatively expensive in computation time. This is expected due to the nonlocal feature of the half-derivative term. We note that this point will be addressed more precisely in the next section.

## 4 | A SECOND NUMERICAL SCHEME

In order to improve the precision and the efficiency of the numerical scheme used before, we construct in this section a second numerical scheme associated to (6) based on a splitting method as described in [26]. In order to construct this scheme, we use the so-called diffusive realization of the half-order derivative. We refer the interested readers to [20, 30, 31]. To this end, we denote by $I^{1/2}u(t) = \frac{1}{\sqrt{\pi}} \int_0^t \frac{u(s)}{\sqrt{t-s}}ds$, the Riemann-Liouville half-order integral and by $D^{1/2}u(t) = \frac{1}{\sqrt{\pi}} \frac{d}{dt} \int_0^t \frac{u(s)}{\sqrt{t-s}}ds$, the Riemann-Liouville half-order derivative. We recall that a diffusive realization of $I^{1/2}u$ is described as follows
<figure><img src="image_3.png" /><figcaption>FIGURE 2 Numerical solutions and the ratios $R_2$ and $R_\infty$ using the Gear scheme when $v = 1$ and $\Delta t = 0.1$ [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 3 Numerical solutions and the ratios $R_2$ and $R_\infty$ using the Gear scheme when $v = 0$ and $\Delta t = 0.1$ [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**FIGURE 4** Numerical solutions and the ratios $R_2$ and $R_\infty$ using the Gear scheme when $\nu = 1$ or 0.1 and $\Delta t = 0.1$
[Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<table><caption>TABLE 2 Decay rates of the solutions when varying the parameters of the Gear scheme with $\Delta t = h = 0.1$</caption><thead><tr><th>Viscosity<br>ν</th><th>Dispersive<br>term β</th><th>Nonlinear<br>term γ</th><th>Diffusive<br>term α</th><th>L<sup>2</sup><br>decay rate</th><th>L<sup>∞</sup><br>decay rate</th><th>Computational<br>time (sec)</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>0</td><td>1</td><td>-0.73</td><td>-0.98</td><td>2214.46</td></tr><tr><td>1</td><td>1</td><td>1</td><td>1</td><td>-0.72</td><td>-0.96</td><td>2215.47</td></tr><tr><td>1</td><td>0</td><td>0</td><td>0</td><td>-0.76</td><td>-1.03</td><td>2258.77</td></tr><tr><td>0</td><td>1</td><td>1</td><td>1</td><td>-0.25</td><td>-0.52</td><td>190.12</td></tr><tr><td>0</td><td>1</td><td>0</td><td>1</td><td>-0.25</td><td>-0.5</td><td>193.49</td></tr><tr><td>0</td><td>0</td><td>1</td><td>1</td><td>-0.25</td><td>-0.5</td><td>191.24</td></tr><tr><td>1</td><td>0.1</td><td>1</td><td>1</td><td>-0.72</td><td>-0.95</td><td>2142.32</td></tr><tr><td>0.1</td><td>1</td><td>0</td><td>1</td><td>-0.79</td><td>-1.12</td><td>2181.37</td></tr><tr><td>1</td><td>0</td><td>0.1</td><td>0.1</td><td>-0.75</td><td>-1.02</td><td>2158.70</td></tr></tbody></table>

$$ \begin{cases} \partial_t \psi(t, \sigma) = -\sigma^2 \psi(t, \sigma) + \frac{2}{\pi} u(t), & \psi(0, \sigma) = 0, \forall \sigma \ge 0, \\ I^{1/2} u(t) = \int_0^{+\infty} \psi(t, \sigma) d\sigma. \end{cases} \tag{28} $$

where $\sigma$ is a new variable not physically relevant. Hence, a diffusive realization of the half-order derivative $D^{1/2}u(t)$ can be deduced by derivation as follows:

$$ \begin{cases} \partial_t \psi(t, \sigma) = -\sigma^2 \psi(t, \sigma) + \frac{2}{\pi} u(t), & \psi(0, \sigma) = 0, \forall \sigma \ge 0, \\ D^{1/2} u(t) = \int_0^{+\infty} \left( \frac{2}{\pi} u(t) - \sigma^2 \psi(t, \sigma) \right) d\sigma. \end{cases} \tag{29} $$

Then we extend the diffusive realization (29) for the functions $u$ depending on time and space as follows.

$$ \begin{cases} \partial_t \psi(t, x, \sigma) = -\sigma^2 \psi(t, x, \sigma) + \frac{2}{\pi} u(t, x), & \psi(0, x, \sigma) = 0, \forall \sigma \ge 0, \\ D^{1/2} u(t, x) = \int_0^{+\infty} \left( \frac{2}{\pi} u(t, x) - \sigma^2 \psi(t, x, \sigma) \right) d\sigma. \end{cases} \tag{30} $$

## 4.1 | Presentation of the model

The nonlocal model (6) can be written as a PDE-ODE coupled system, using (30), as follows

$$ \left\{ \begin{aligned} \partial_t u(t,x) + \partial_x(u + \frac{\gamma}{2}u^2) &= -\sqrt{\nu} \int_0^\infty \left( \frac{2}{\pi}u(t,x) - \sigma^2 \psi(t,x,\sigma) \right) d\sigma \\ &\quad + \alpha u_{xx}(t,x) + \beta u_{txx}(t,x), \quad t > 0, x \in \mathbb{R} \\ \partial_t \psi(t,x,\sigma) &= -\sigma^2 \psi(t,x,\sigma) + \frac{2}{\pi}u(t,x), \quad t > 0, x \in \mathbb{R}, \sigma \ge 0, \end{aligned} \right. \tag{31} $$

supplemented with the initial conditions

$$ \begin{cases} \forall x \in \mathbb{R}, \forall \sigma \ge 0, \quad \psi(0, x, \sigma) = 0, \\ \forall x \in \mathbb{R}, \quad u(0, x) = u_0(x). \end{cases} $$

In order to approximate the Riemann-Liouville half-order derivative in (6), we need to approximate the generalized integral in (31). To this end, we use a quadrature formula with $N_m$ points. We note
by $w_i$ the weights and by $\sigma_i$ the nodes (or abscissae) of the appropriate quadrature method used in the approximation. We get

$$ 
\begin{aligned}
D^{1/2} u(t, x) &\approx \sum_{i=1}^{N_m} w_i \left( \frac{2}{\pi} u(t, x) - \sigma_i^2 \psi(t, x, \sigma_i) \right) \\
&= \sum_{i=1}^{N_m} w_i \left( \frac{2}{\pi} u(t, x) - \sigma_i^2 \psi_i(t, x) \right),
\end{aligned}
$$

Hence, the system (31) is written as a first order system as follows

$$ 
\left\{
\begin{array}{l}
\displaystyle \partial_t u(t, x) + \partial_x \left(u + \frac{\gamma}{2} u^2\right) = -\sqrt{\nu} \sum_{i=1}^{N_m} w_i \left(\frac{2}{\pi} u(t, x) - \sigma_i^2 \psi_i(t, x)\right) \\
\displaystyle \qquad + \alpha u_{xx}(t, x) + \beta u_{txx}(t, x), \quad t > 0, x \in \mathbb{R}, \\
\\
\forall i = 1, \dots, N_m, \\
\\
\displaystyle \partial_i \psi_i(t, x) = -\sigma_i^2 \psi_i(t, x) + \frac{2}{\pi} u(t, x), \quad t > 0, x \in \mathbb{R},
\end{array}
\right.
\quad (32)
$$

endowed with the initial conditions

$$ 
\begin{array}{l}
\forall x \in \mathbb{R}, \forall i = 1, \dots, N_m, \quad \psi_i(0, x) = 0, \\
\forall x \in \mathbb{R}, \quad u(0, x) = u_0(x).
\end{array}
$$

Now, we note by

$$ U = (u, \psi_1, \dots, \psi_{N_m})^T, $$

the vector of $(N_m + 1)$ unknowns, by

$$ \mathcal{F}(U) = (u + \frac{\gamma}{2} u^2, 0, \dots, 0)^T, $$

and finally

$$ 
S = \begin{pmatrix}
-\sqrt{\nu} \sum_{i=1}^{N_m} w_i & \sqrt{\nu} w_1 \sigma_1^2 & \cdots & \cdots & \sqrt{\nu} w_{N_m} \sigma_{N_m}^2 \\
\frac{2}{\pi} & -\sigma_1^2 & 0 & \cdots & 0 \\
\frac{2}{\pi} & 0 & -\sigma_2^2 & \cdots & 0 \\
\vdots & \vdots & \ddots & \ddots & \vdots \\
\frac{2}{\pi} & 0 & \cdots & \cdots & -\sigma_{N_m}^2
\end{pmatrix}.
$$

Thus, the problem (32) is written in the following form

$$ \partial_t U + \partial_x \mathcal{F}(U) = S(U) + G_1 \partial_x^2 U - G_2 \partial_t \partial_x^2 U, \quad (33) $$

where $G_1$ and $G_2$ are diagonal matrices of order $N_m + 1$. In the sequel, we introduce the so-called splitting scheme.
## 4.2 | The splitting method

Let $\Delta t > 0$, for all $n \ge 0$, we recall that $t^n = n\Delta t$ and

$$U^n(x) \approx U(n\Delta t, x).$$

From Equation 33, we consider the propagation equation

$$\partial_t U + \partial_x \mathcal{F}(U) = \mathbf{G}_1 \partial_x^2 U - \mathbf{G}_2 \partial_t \partial_x^2 U, \quad (34)$$

and the diffusive equation

$$\partial_t U = \mathbf{S}(U). \quad (35)$$

We note by $\mathbf{H}_a$ (respectively, $\mathbf{H}_b$) the discrete operator of the solution of (34) (respectively, the solution of (35)). Then a Strang Splitting method of order 2 ([32, 33]) between $t_n$ and $t_{n+1}$ is used to solve respectively (34) and (35) as follows

$$\begin{align}
U^{(1)} &= \mathbf{H}_b \left( \frac{\Delta t}{2} \right) U^n, \notag \\
U^{(2)} &= \mathbf{H}_a (\Delta t) U^{(1)}, \notag \\
U^{n+1} &= \mathbf{H}_b \left( \frac{\Delta t}{2} \right) U^{(2)}. \tag{36}
\end{align}$$

Here, the constructed operators $\mathbf{H}_a$ and $\mathbf{H}_b$ are stable and of order 2. Then, the scheme (36) provides an approximation of order 2 in time to the problem (33).

In the sequel, we present the discretization of (34) and (35).

The propagation equation (34). Here $u$ is a solution of the BBM equation. We perform a semidiscrete in time scheme: we use a Crank-Nicolson scheme for the linear part and Adams-Bashforth scheme (see [27]) for the nonlinear part. For the space discretization, we use standard Fourier methods.

First, we note that the first approximate solution $\hat{u}^1$ is performed using a fixed-point method that verifies the semidiscrete scheme (of order 2).

$$(1 + \beta \xi^2) \frac{\hat{u}^1 - \hat{u}^0}{\Delta t} = \frac{\hat{u}^1 + \hat{u}^0}{2} (-\alpha \xi^2 - i\xi) - \frac{i\gamma\xi}{8} \left( (\hat{u}^0)^2 + (\hat{u}^1)^2 \right). \quad (37)$$

Then for $n \ge 1$, the discrete scheme is given by

$$(1 + \beta \xi^2) \frac{\hat{u}^{n+1} - \hat{u}^n}{\Delta t} = \frac{\hat{u}^{n+1} + \hat{u}^n}{2} (-\alpha \xi^2 - i\xi) - \frac{i\gamma\xi}{8} (3(\hat{u}^n)^2 - (\hat{u}^{n-1})^2). \quad (38)$$

Diffusion equation (35). We can solve mathematically (35) as follows

$$\mathbf{H}_b \left( \frac{\Delta t}{2} \right) U = e^{s \frac{\Delta t}{2}} U. \quad (39)$$
In addition, for all $N_m > 0$, the exponential of the matrix $S^{\Delta t/2}$ is calculated numerically using the “scaling and squaring” method with a (6, 6) Padé approximation [34] (corresponding with Matlab® to the function “expm(S^{\Delta t/2}”).). We recall that the $(p, q)$ Padé approximation of $e^A$ is given by

$$R_{pq}(A) = (D_{pq}(A))^{-1}N_{pq}(A),$$

where

$$N_{pq}(A) = \sum_{j=0}^{p} \frac{(p+q-j)!p!}{(p+q)!j!(p-j)!} A^j,$$

$$D_{pq}(A) = \sum_{j=0}^{p} \frac{(p+q-j)!q!}{(p+q)!j!(q-j)!} (-A)^j.$$

Hence, the numerical scheme (36) is a three-step scheme.

We note that the scheme (38) as well as the Strang Splitting method are of order two. It follows that (36) is of order 2. This will be numerically verified in the sequel.

## 4.3 | Quadrature method

In this subsection, we are interested in the choice of the $2N_m$ coefficients $w_i$ and $\sigma_i$ of the diffusive representation given in (32). These coefficients aim to approach the improper integrals in the form

$$\int_{0}^{+\infty} \psi(\sigma) d\sigma \approx \sum_{i=1}^{N_m} w_i \psi(\sigma_i). \qquad (40)$$

For seek of convenience, we drop here the variables $t$ and $x$. The choice of these coefficients is an important issue, because it affects directly the accuracy of the method and the efficiency of the approximation. We note that many methods developed in the literature are especially based on the Gauss quadrature. The choice of the quadrature method was thoroughly discussed in [30] when considering the KdV equation with the Riemann-Liouville half-derivative. Based on this work, we present the quadrature method that will be used in the remaining of this article.

**Gauss-Jacobi quadrature method.** Here, we aim to approximate the generalized integral (40) with the Gauss-Jacobi quadrature method. It consists in transforming the domain of integration from $[0, +\infty[$ to $[-1, 1]$ then applying the Gauss-Jacobi quadrature to the resulting integral. Hence, we choose the following change of variables

$$\sigma = \frac{1-z}{1+z} \quad \text{then} \quad \frac{d\sigma}{dz} = \frac{-2}{(1+z)^2}.$$

It follows that

$$\int_{0}^{+\infty} \psi(\sigma) d\sigma = \int_{-1}^{1} \frac{2}{(1+z)^2} \psi\left(\frac{1-z}{1+z}\right) dz.$$

Moreover, letting

$$\tilde{\psi}(z) = \frac{2}{(1+z)^2} \psi\left(\frac{1-z}{1+z}\right),$$
<table><caption>**TABLE 3** Weights and abscissae of the Gauss-Laguerre and Gauss-Jacobi quadrature formulas for $N_m = 8$</caption><thead><tr><th colspan="3">Gauss-Jacobi</th></tr><tr><th>w<sub>i</sub></th><th>z<sub>i</sub></th><th>σ<sub>i</sub></th></tr></thead><tbody><tr><td>128.3897</td><td>-0.9603</td><td>49.3650</td></tr><tr><td>10.7575</td><td>-0.7967</td><td>8.8361</td></tr><tr><td>2.7870</td><td>-0.5255</td><td>3.2153</td></tr><tr><td>1.0879</td><td>-0.1834</td><td>1.4493</td></tr><tr><td>0.5179</td><td>0.1834</td><td>0.6899</td></tr><tr><td>0.2696</td><td>0.5255</td><td>0.3110</td></tr><tr><td>0.1378</td><td>0.7967</td><td>0.1132</td></tr><tr><td>0.0527</td><td>0.9603</td><td>0.0203</td></tr></tbody></table>

we get

$$ \int_{0}^{+\infty} \psi(\sigma) d\sigma = \int_{-1}^{1} \tilde{\psi}(z) dz \approx \sum_{i=0}^{N_m} \mu_i \tilde{\psi}(z_i), $$

where $\mu_i$ (resp. $z_i$) are the weights (resp. the nodes) of the standard Gauss-Jacobi quadrature formula over [-1, 1]. For illustrative purposes, we present in Table 3 the weights and the nodes of this quadrature formula for $N_m = 8$. By identification, we deduce that the quadrature coefficients in (40) are given by

$$ w_i = \frac{2}{(1 + z_i)^2} \mu_i, \quad \sigma_i = \frac{1 - z_i}{1 + z_i}. $$

We note that this strategy was proposed by Diethelm in [8] for the approximation of Caputo fractional-order derivative.

## 4.4 | Numerical results and discussion.

We begin with justifying the convergence in time of the splitting scheme (37)–(38). To this end, we choose $h = 0.1$, $L = 800$, $N_m = 20$, $T = 100$, and $\alpha = \beta = \gamma = \nu = 1$. We denote by $u_{\text{Ref}}^n$ the reference solution when $\Delta t = 0.05$ and $u^n$ the numerical solution for different time steps $\Delta t$. Also, we denote by $E^n(\Delta t) = \|u_{\text{Ref}}^n - u^n\|_2$ the $L^2$-error in terms of $\Delta t$. The results are presented in Figure 5. We see that the error $E^n$ decreases when the time step $\Delta t$ decreases. Also, we may determine numerically the order of the scheme (37)–(38). In fact, the measured values are close to a straight line with slope 2 which means that $E^n(\Delta t) \approx C \Delta t^2$ where $C$ is a constant. We conclude that the numerical scheme (37)–(38) is of order 2 in time.

Then, we aim to study numerically the convergence of the splitting scheme (37)–(38) with respect to the number of Gauss-Jacobi quadrature points. To this end, we consider the parameters $\Delta t = h = 0.1$, $L = 800$, $T = 100$ and $\nu = \alpha = \beta = \gamma = 1$. We denote by $u_{\text{Ref}}^n$ the reference solution for $N_m = 20$ and by $u^n$ the numerical solution for different values of $N_m$. Some other numerical tests show that convergence is obtained with $N_m = 20$ (see also [26]). Moreover, we denote by $E^n(N_m) = \|u_{\text{Ref}}^n - u^n\|_2$ the $L^2$-error in terms of quadrature points. We plot, in Figure 6, the Error $E^n$ with respect to $N_m$. We see that the measured values are close to a straight line with slope -11, which means that $E^n(N_m) \approx C(\frac{1}{N_m})^{11}$ where $C$ is a constant.
<figure><img src="image_6.png" /><figcaption>**FIGURE 5** Error of the time discretization using the splitting scheme [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

Besides, since the error is less than $10^{-5}$ from $N_m = 15$, we choose this number of points to realize the simulations in the sequel.

Now, we examine the convergence of the decay rates of the Gear scheme (25)–(26) and the splitting scheme (37)–(38). To this end, we take the parameters values $L = 800$, $\Delta t = h = 0.1$, $\alpha = \beta = 1$, $\nu = 0.1$, and $\gamma = 0.5$. We calculate the solutions up to time $T = 100$. We denote by $R_p^{Gr}$ the ratio in norm $L^p$ of the Gear scheme (25)–(26) and by $R_p^{Sp}$ the ratio in norm $L^p$ of the splitting scheme (37)–(38) where $p = 2$ or $\infty$. Also, we denote by $\text{Err}_p = R_p^{Sp} - R_p^{Gr}$ the difference with respect to the time of the decay rate obtained by the two approximations. In Figure 7, we plot $\text{Err}_2$ and $\text{Err}_\infty$ with respect to time $t$. We observe that both errors $\text{Err}_2$ and $\text{Err}_\infty$ are decreasing to zero especially for large times for all the time step $\Delta t$ chosen, that is the main part of interest of our study. We conclude that the numerical schemes converges for large time.

<figure><img src="image_7.png" /><figcaption>**FIGURE 6** Error in terms of $N_m$ using the Splitting scheme [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**FIGURE 7** Differences of the decay rates between the splitting and the Gear schemes when $N_m = 15$, $\alpha = \beta = 1$, $\nu = 0.1$, $\gamma = 0.5$, and $\Delta t = 0.1$ [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

In the sequel, we study the time convergence of the decay rates $R_2$ and $R_\infty$ of the splitting scheme (37)–(38). To this end, we fix the parameters values to $h = 0.1$, $L = 800$, $T = 100$, $\alpha = 2$, $\beta = \gamma = 1$, $\nu = 0.5$. We denote by $R_p^{\text{Ref}}$ the ratio in norm $L^p$ when the time step $\Delta t = 0.05$ and by $R_p^n$ the ratio in norm $L^p$ for different values of $\Delta t$ when $p = 2$ or $\infty$. Also, we denote by $\text{Err}_p^n = R_p^{\text{Ref}} - R_p^n$ the error in time of the decay rates. In Figure 8, we present the errors $\text{Err}_2^n$ and $\text{Err}_\infty^n$ when increasing the time step $\Delta t$ from 0.05 to 0.25. We observe that both errors decreases when increasing the time step $\Delta t$. Moreover, we observe that the errors are less than $10^{-3}$ for large times when $\Delta t = 0.1$. Hence, we conclude that this time step is suffisant to get accurate numerical results.

Finally, we determine the computation time of solutions (in seconds) using the splitting scheme (37)–(38) for different parameters values when $\Delta t = h = 0.1$. Results are displayed in Table 4. As it is expected, the time elapsed to calculate the numerical solution using the splitting scheme is reduced. It represents 1/5 of the time necessary to calculate solutions with the Gear scheme (25)–(26). We deduce that the use of the splitting scheme performs a numerical solution more accurate and in a relatively shorter time of computation.
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**FIGURE 8** Differences of the decay rates between a solution of reference with $\Delta t = 0.05s$ and solutions obtained with larger time steps, using splitting scheme where $N_m = 15$, $\alpha = 2$, $\nu = 0.5$, and $\gamma = \beta = 1$ [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<table><caption>**TABLE 4** Decay rates of the solutions when varying parameters of the Splitting scheme with $\Delta t = h = 0.1$</caption><thead><tr><th>Viscosity<br>ν</th><th>Dispersive<br>term β</th><th>Nonlinear<br>term γ</th><th>Diffusive<br>term α</th><th>L<sup>2</sup><br>decay rate</th><th>L<sup>∞</sup><br>decay rate</th><th>Computational<br>time (sec)</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>0</td><td>1</td><td>-0.74</td><td>-0.98</td><td>400.95</td></tr><tr><td>1</td><td>1</td><td>1</td><td>1</td><td>-0.73</td><td>-0.97</td><td>401.8</td></tr><tr><td>1</td><td>0</td><td>0</td><td>0</td><td>-0.76</td><td>-1.03</td><td>408.67</td></tr><tr><td>0</td><td>1</td><td>1</td><td>1</td><td>-0.25</td><td>-0.52</td><td>190.12</td></tr><tr><td>0</td><td>1</td><td>0</td><td>1</td><td>-0.25</td><td>-0.5</td><td>193.49</td></tr><tr><td>0</td><td>0</td><td>1</td><td>1</td><td>-0.25</td><td>-0.5</td><td>191.24</td></tr><tr><td>1</td><td>0.1</td><td>1</td><td>1</td><td>-0.73</td><td>-0.96</td><td>399.78</td></tr><tr><td>0.1</td><td>1</td><td>0</td><td>1</td><td>-0.84</td><td>-1.16</td><td>394.59</td></tr><tr><td>1</td><td>0</td><td>0.1</td><td>0.1</td><td>-0.76</td><td>-1.02</td><td>393.23</td></tr></tbody></table>
# 5 | CONCLUSION

In this article, we have constructed two numerical schemes to approximate the solutions and the decay rates to an asymptotical water wave model where the nonlocal viscous term is described by the Riemann-Liouville half derivative. We compare our numerical results to those given in [1, 2, 10]. We show that using the diffusive realization of the nonlocal operator supplemented with a splitting scheme leads to a very interesting gain of time computing of the numerical solution. This numerical scheme enables us to approximate nonlocal models in a shorter time when comparing with classical methods as the Gear scheme. A challenging issue is to address analytically the asymptotical behavior of the initial value problem. This question will be the subject of a future work.

## ACKNOWLEDGMENTS

The authors would like to express their appreciation to Referees for their constructive remarks and suggestions. Also, the authors are grateful to Professor Olivier Goubet for his useful comments on section 2. This work was partially supported by a grant from the *Simons Foundation*.

## ORCID

Imen Manoubi http://orcid.org/0000-0003-2243-1330

## REFERENCES

[1] S. Dumont, J.-B Duval, *Numerical investigation of asymptotical properties of solutions to models for waterwaves with nonlocal viscosity*, Int. J. Numer. Anal. Model. vol. 10 (2013) pp. 333–349.

[2] I. Manoubi, *Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative*, Discrete Continuous Dyn. Syst. vol. 19 (2014) pp. 2837–2863.

[3] J. L. Bona, M. Chen, J.-C Saut, *Boussinesq equations and other systems for small-amplitude long waves in nonlinear dispersive media I: Derivation and the linear theory*, J. Nonlinear Sci. vol. 12 (2002) pp. 283–318.

[4] T. Kakutani, M. Matsuuchi, *Effect of viscosity on long gravity waves*, J. Phys. Soc. Jpn. vol. 39 (1975) pp. 237–246.

[5] E. Ott, R. N. Sudan, *Damping of solitary waves*. Phys. Fluids vol. 13 (1970) pp. 1432–1434.

[6] P. Liu, A. Orfila, *Viscous effects on transient long wave propagation*, J. Fluid Mech. vol. 520 (2004) pp. 83–92.

[7] D. Dutykh, F. Dias, *Viscous potential free surface flows in a fluid layer of finite depth*, C.R.A.S, Série I vol. 345 (2007) pp. 113–118.

[8] K. Diethelm, *An investigation of some nonclassical methods for the numerical approximation of Caputo-type fractional derivatives*, Numer. Algorithms vol. 47 (2008) pp. 361–390.

[9] M. Chen, S. Dumont, L. Dupaigne, O. Goubet, *Decay of solutions to a water wave model with a nonlocal viscous dispersive term*, Discrete Continuous Dyn. Syst. vol. 27 (2010) pp. 1473–1492.

[10] O. Goubet, I. Manoubi, *Theoretical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach*, Adv. Nonlinear Anal. (2017). https://doi.org/10.1515/anona-2016-0274.

[11] L. Di Menza, *Approximations numériques d'équations de Schrödinger non linéaires et de modèles associés*, PhD Thesis, Bordeaux 1, France, 1995.

[12] A. C. Galucio, J.-F Deü, F. Dubois, *The $G^a$-scheme for approximation of fractional derivatives: application to the dynamics of dissipative systems*, J. Vib. Control vol. 14 (2008) pp. 1597–1605.

[13] A. C. Galucio, J.-F Deü, S. Mengu, F. Dubois, *An adaptation of the Gear scheme for fractional derivatives*, Comput. Methods Appl. Mech. Eng. vol. 195 (2006) pp. 6073–6085.

[14] C. Koh, J. Kelly, *Application of fractional derivatives to seismic analysis of base-isolated models*, Earthq. Eng. Struct. Dyn. vol. 19 (1990) pp. 229–241.
[15] J.-R Li, *A fast time stepping method for evaluating fractional integrals*, SIAM J. Sci. Comput. vol. 31 (2010) pp. 4696-4714.

[16] A. Shokooh, L. Suarez, *A comparison of numerical methods applied to a fractional model of damping materials*, J. Vib. Control vol. 5 (1999) pp. 331-354.

[17] W. Zhang, N. Shimizu, *Numerical algorithm for dynamic problems involving fractional operators*, JSME Int. J. Ser. C vol. 41 (1998) pp. 364-370.

[18] G. Montseny, *Diffusion monodimensionnelle et intégration d'ordre 1/2*. *Internal LAAS Report N.* 91232, 1991.

[19] G. Montseny, J. Audounet, B. Mbodge, *Modèle simple d'amortisseur viscoélastique. Application à une corde vibrante, in Lecture notes in Control and Information Sciences* Eds. R. F. Curtain, A. Bensoussan, J. L. Lions-Springer Verlag 1993;185:436-446.

[20] G. Montseny, J. Audounet, B. Mbodge, *Optimal models of fractional integrators and application to systems with fading memory*, *IEEE International Conference on Systems, Man and Cybernetics, Le Touquet France* 1993; 17-20:65-70.

[21] O. Staffans, *Well-posedness and stabilizability of a viscoelastic equation in energy space*, Trans. Am. Math. Soc. vol. 345 (1994) pp. 527-575.

[22] J. Audounet, V. Giovangigli, J. Roquejoffre, *A threshold phenomenon in the propagation of a point-source initiated flame*, Physica D vol. 121 (1998) pp. 295-316.

[23] T. Helie, D. Matignon, *Diffusive representations for the analysis and simulation of flared acoustic pipes with visco-thermal loss*, Math. Models Methods Appl. Sci. vol. 16 (2006) pp. 503-536.

[24] D. Matignon, C. Prieur, *Asymptotic stability of linear conservative systems when coupled with diffusive systems*, ESAIM: Control Optim. Calc. Var. vol. 11 (2005) pp. 487-507.

[25] G. Montseny, J. Audounet, D. Matignon. *Diffusive representation for pseudodifferentially damped nonlinear systems*, Nonlinear Control in the year 2000 2000; 2:163-182.

[26] B. Lombard, J-F Mercier, *Numerical modeling of nonlinear acoustic waves in a tube connected with Helmholtz*, J. Comput, Phys. vol. 259 (2014) pp. 421-443.

[27] H. Kalisch, J. Bona, *Models for internal waves in deep water*, Discrete Continuous Dyn. Syst. vol. 6 (2000) pp. 1-20.

[28] D. Dutykh, *Viscous-potential free-surface flows and long wave modelling*, Eur. J. Mech. B Fluids vol. 28 (2009) pp. 430-443.

[29] C. J. Amick, J. L. Bona, M. E. Schonbek, *Decay of solutions of some nonlinear wave equations*, J. Differ. Equ. vol. 81 (1989) pp. 1-49.

[30] I. Manoubi, *Modèle visqueux asymptotique pour la propagation d'une onde dans un canal*, PhD Thesis, Monastir, Tunisia, 2014.

[31] G. Montseny, *Diffusive representation of pseudo-differential time operators*, ESAIM Proc Fract. Differ. Syst. Models Methods Appl. vol. 5 (1998) pp. 159-175.

[32] H. Holden, K. Karlsen, N. Risebro, T. Tao, *Operator splitting for the KDV equation*, Math. Comput. vol. 80 (2011) pp. 821-846.

[33] R. LeVeque, *Numerical methods for conservation laws*, 2nd ed., Birkhäuser-Verlag, Basel, Boston, Berlin, 1992.

[34] C. Moler, C. Van Loan, *Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later*, SIAM Rev. vol. 45 (2003) pp. 3-49.