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

HAL
open science

Numerical analysis of a water wave model with a
nonlocal viscous dispersive term using the diffusive
approach

Serge Dumont, Imen Manoubi

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

To cite this version:

Serge Dumont, Imen Manoubi. Numerical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach. Mathematical Methods in the Applied Sciences, 2018, 41 (12), pp.4810-4826. 10.1002/mma.4932 . hal-01818120

HAL Id: hal-01818120

https://hal.science/hal-01818120

Submitted on 17 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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach

S. Dumont<fnref n="1" />,<fnref n="2" /> | I. 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

PHC Utique ASEO

MSC Classification: 35Q35; 35Q53

In this paper, we numerically study the water wave model with a nonlocal viscous term

$$u_t + u_x + \beta u_{xxx} + \sqrt{v} D^{1/2} u(t) + \gamma u u_x = \alpha u_{xx},$$

where $D^{1/2}u(t) = \frac{1}{\sqrt{\pi}} \frac{\partial}{\partial t} \int_0^t \frac{u(s)}{\sqrt{t-s}} ds$ is the Riemann-Liouville half-order derivative in time. We propose and compare different numerical schemes based on the diffusive realizations of fractional operators.

## KEYWORDS

decay rate, diffusive realization, fractional derivatives, nonlocal viscous model, time splitting, water waves

# 1 | INTRODUCTION

## 1.1 | State of the art

Modeling the effects of viscosity in the propagation of long waves is a challenging issue. In the last decade, Liu-Orfila<fnref n="1" /> and Dutykh-Dias<fnref n="2" /> have independently derived mathematical models in 2D and in 3D where the viscosity is modeled with nonlocal in time operators. In the case of a mono-dimension wave, the model is reduced as follows:

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

where $\frac{1}{\sqrt{\pi}} \int_0^t \frac{u_t(s)}{\sqrt{t-s}} ds$ represents the Caputo half derivative in time. Here, $u$ is the horizontal velocity of the fluid, $-\alpha u_{xx}$ is the usual diffusion, $\beta u_{xxx}$ is the geometric dispersion, and $\frac{1}{\sqrt{\pi}} \int_0^t \frac{u_t(s)}{\sqrt{t-s}} ds$ stands for the nonlocal diffusive-dispersive term. We denote as $\beta$, $v$, and $\alpha$ the parameters dedicated to balance or unbalance the effects of viscosity and dispersion against nonlinear effects. Some numerical results of the decay rate of a solution of (1) are performed in Chen et al<fnref n="3" /> and Dumont and Duval.⁴ Also, in the recent work,⁵ the author has considered the model (1) where she used the Riemann-Liouville half derivative instead of that of Caputo, namely,

$$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 + uu_x = \alpha u_{xx}. \quad (2)$$

<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>
Comparing with the problem (1), she improved the estimates on the decay rates of the solution. Particularly, she proved that the local and global existence result and decay estimates for the integro-differential equation 2 when $\beta = 0$, $\nu = \alpha = 1$ supplemented with the initial condition $u_0 \in L^1(\mathbb{R}) \cap L^2(\mathbb{R})$. Precisely, she stated the following theorem.

**Theorem 1.** 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 \otimes 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 $\otimes$ is the time-space convolution product defined by

$$v \otimes w(t, x) = \int_{0}^{t} \int_{\mathbb{R}} v(t-s, x-y)w(s, y)dsdy$$

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 Manoubi.⁵ The approximation of time-fractional operators has received a lot of interest during last decades for their wide application in fluid and solid mechanics and visco-elasticity. The formulation of a numerical stable scheme is crucial but also a difficult question because of the nonlocal feature of such operators. The classical methods used in the literature⁶⁻¹¹ consist on the approximation of these fractional operators using either convolution integrals or the so-called Gear scheme for fractional operators. However, these methods are expensive in terms of run time evaluations and of the memory required to store the function values, especially for large time simulations. Also, Zang and Xu proposed in Zhang and Xu¹² numerical schemes based on finite difference/spectral approximations to a water wave model. In particular, they used the known $(2 - \alpha)$ order scheme for the $\alpha$-order fractional derivative and a mixed linearization for the nonlinear term. They proved the unconditionally stability and some convergence rates for the schemes. Moreover, in their work,¹³ Liu et al constructed finite difference/finite element method for a nonlinear time-fractional fourth-order reaction-diffusion problem. They proved unconditionally stability of the schemes and some a priori estimates of $L^2$ norm with optimal order of convergence. Independently, a number of authors suggested a method for the direct approximation of the fractional operators. This method, called diffusive representation,¹⁴ is an operator theory approach and is devoted to *sech* pseudo-differential equations and time nonlocal problems. This approach was developed by Montseny,¹⁵ Montseny et al,¹⁶,¹⁷ and Staffans.¹⁸ Different applications of this approach can be found in other studies.¹⁹⁻²² 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, which has many advantages from the mathematical and numerical points of view.

## 1.2 | Diffusive formulation of the model

In the literature, there are several diffusive realizations of the Riemann-Liouville half-order derivative $D^{1/2}u(t)$. We recall in the following some of these formulations. First, the diagonal form of the diffusive realization of $D^{1/2}u(t)$ is given for all $t > 0$ by

$$\left\{
\begin{array}{l}
\partial_t \psi(t, \sigma) = -\sigma \psi(t, \sigma) + u(t), \quad \psi(0, \sigma) = 0, \quad \sigma \ge 0, \\
D^{1/2} u(t) = \int_0^{+\infty} \frac{1}{\pi \sqrt{\sigma}} \partial_t \psi(t, \sigma) d\sigma.
\end{array}
\right.
\quad (5)$$
Second, the partial differential equation (PDE) form of the diffusive realization of $D^{1/2}u(t)$ is given for all $t > 0$ by

$$ \left\{ \begin{aligned} \partial_t \Phi(t, y) &= \Phi_{yy}(t, y) + u(t) \otimes \delta_{y=0}, \quad \Phi(0, y) = 0, \quad y \in \mathbb{R}, \\ D^{1/2} u(t) &= 2 < \delta_{y=0}, \partial_t \Phi(t, y) > D', \quad D' = 2 \frac{d}{dt} \Phi(t, 0), \end{aligned} \right. \qquad (6) $$

where $\delta_{y=0}$ is the Dirac delta function at $y = 0$ and $u(t) \otimes \delta_{y=0}$ is the tensorial product in the distributions sense of the applications $t \mapsto u(t)$ and $y \mapsto \delta_{y=0}$. Finally the diffusive realization of $D^{1/2}u(t)$ used in the remaining of this article is given for all $t > 0$ by

$$ \left\{ \begin{aligned} \partial_t \phi(t, \sigma) &= -\sigma^2 \phi(t, \sigma) + \frac{2}{\pi} u(t), \quad \phi(0, \sigma) = 0, \quad \sigma \ge 0, \\ D^{1/2} u(t) &= \int_0^{+\infty} \left( \frac{2}{\pi} u(t) - \sigma^2 \phi(t, \sigma) \right) d\sigma. \end{aligned} \right. \qquad (7) $$

We note that the author has used the diffusive realizations (6) and (7) in her PhD thesis to study mathematically the integro-differential equation 2 when $\beta = 0, \nu = \alpha = 1$ using 2 different methods. For more details, we refer the readers to Manoubi.²⁵ In this paper, we study numerically the integro-differential equation 2 using the diffusive approach. Hence, we extend the diffusive realization (7) of the Riemann-Liouville half-order derivative as follows:

$$ \left\{ \begin{aligned} \partial_t \psi(t, x, \sigma) &= -\sigma^2 \psi(t, x, \sigma) + \frac{2}{\pi} u(t, x), \quad \psi(0, x, \sigma) = 0, \quad \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{aligned} \right. \qquad (8) $$

Then, the integro-differential equation 2 can be written as a PDE-ordinary differential equation (ODE) coupled system as follows:

$$ \left\{ \begin{aligned} u_t(t, x) + \sqrt{\nu} \int_0^{+\infty} \left( \frac{2}{\pi} u(t, x) - \sigma^2 \psi(t, x, \sigma) \right) d\sigma = \alpha u_{xx}(t, x) \\ -u_x(t, x) - \beta u_{xxx}(t, x) - \gamma u(t, x)u_x(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. \qquad (9) $$

supplemented with the initial conditions

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

The remaining of this article is as follows. Different quadrature methods are presented in Section 2 to approximate the Riemann-Liouville half-order derivative. In Section 3, we present a first numerical scheme associated to the system (9) where we used the quadrature methods followed by several numerical simulations and discussion. In Section 4, we present a second numerical scheme for (9) using a splitting procedure. Numerical results are proposed in order to check the efficiency of the numerical scheme. A dispersion analysis is presented in Section 5.

## 2 | QUADRATURE METHODS

To approximate the Riemann-Liouville half-order derivative in (2), we need to approximate the generalized integral in (9). To this end, we use a quadrature formula with $N$ 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} w_i \left( \frac{2}{\pi} u(t, x) - \sigma_i^2 \psi(t, x, \sigma_i) \right) \\ &= \sum_{i=1}^{N} w_i \left( \frac{2}{\pi} u(t, x) - \sigma_i^2 \psi_i(t, x) \right), \end{aligned} $$
where $\psi_i$ verifies the ODE (8). Hence, system (9) is approximated as follows:

$$ egin{cases} u_t(t,x) + \sqrt{v} \sum_{i=1}^{N} w_i \left( \frac{2}{\pi} u(t,x) - \sigma_i^2 \psi_i(t,x) \right) = \alpha u_{xx}(t,x) \\ -u_x(t,x) - \beta u_{xxx}(t,x) - \gamma u(t,x) u_x(t,x), \quad t > 0, x \in \mathbb{R}, \\ \partial_t \psi_i(t,x) = -\sigma_i^2 \psi_i(t,x) + \frac{2}{\pi} u(t,x), \quad t > 0, x \in \mathbb{R}, i = 0, \dots, N, \end{cases} \quad (10) $$

supplemented with the initial conditions

$$ \begin{gathered} \forall x \in \mathbb{R}, \forall i = 0, \dots, N \qquad \psi_i(0,x) = 0, \\ \forall x \in \mathbb{R}, u(0,x) = u_0(x). \end{gathered} $$

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

$$ \int_{0}^{+\infty} \psi(\sigma) d\sigma \approx \sum_{i=1}^{N} w_i \psi(\sigma_i). \quad (11) $$

For seek of convenience, we drop here the variables $t$ and $x$. The choice of the coefficients $w_i$ and $\sigma_i$ is an important issue, because it effects 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. In the sequel, we present some quadrature methods that will be used in the remaining of this article.

## 2.1 | The Gauss-Laguerre quadrature

A first method consists in approximating the generalized integral (11) using the standard Gauss-Laguerre quadrature over the interval $[0, +\infty[$ with the weight function $w(\sigma) = e^{-\sigma}$ (see Abramowitz and Stegun$^{23}$)

$$ \int_{0}^{+\infty} e^{-\sigma} \psi(\sigma) d\sigma \approx \sum_{i=0}^{N} \mu_i \psi(z_i). \quad (12) $$

Particularly, this method is more suitable in the case of an integrand with an exponential decay. The coefficients $\mu_i$ (resp $z_i$) are the weights (resp the nodes) of the Gauss-Laguerre quadrature formula. Hence, (12) implies

$$ \int_{0}^{+\infty} \psi(\sigma) d\sigma \approx \sum_{i=0}^{N} \mu_i e^{z_i} \psi(z_i), $$

As a consequence, the coefficients in (11) are given by

$$ w_i = \mu_i e^{z_i}, \quad \sigma_i = z_i. $$

For illustrative purposes, we present in Table 1 the weights and the nodes of this quadrature method when $N = 8$. We note that this quadrature method was proposed by Yuan and Agrawal$^{24}$ in order to approximate the Caputo fractional derivative of order $\alpha$ where $0 < \alpha < 1$.

## 2.2 | The Gauss-Jacobi quadrature method

In this paragraph, a second method is proposed, approximating the generalized integral (11) with the Gauss-Jacobi quadrature method. For that purpose, it is necessary to transform the domain of integration from $[0, +\infty[$ to $[-1, 1]$, then to apply 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. $$
<table><caption>**TABLE 1** Weights and abscissae of the Gauss-Laguerre and Gauss-Jacobi quadrature formulas for $N = 8$</caption><thead><tr><th colspan="2">Gauss-Laguerre</th><th colspan="3">Gauss-Jacobi</th></tr><tr><th>w<sub>i</sub></th><th>σ<sub>i</sub></th><th>w<sub>i</sub></th><th>z<sub>i</sub></th><th>σ<sub>i</sub></th></tr></thead><tbody><tr><td>0.4377</td><td>0.1703</td><td>128.3897</td><td>-0.9603</td><td>49.3650</td></tr><tr><td>1.0339</td><td>0.9037</td><td>10.7575</td><td>-0.7967</td><td>8.8361</td></tr><tr><td>1.6697</td><td>2.2511</td><td>2.7870</td><td>-0.5255</td><td>3.2153</td></tr><tr><td>2.3769</td><td>4.2667</td><td>1.0879</td><td>-0.1834</td><td>1.4493</td></tr><tr><td>3.2085</td><td>7.0459</td><td>0.5179</td><td>0.1834</td><td>0.6899</td></tr><tr><td>4.2686</td><td>10.7585</td><td>0.2696</td><td>0.5255</td><td>0.3110</td></tr><tr><td>5.8181</td><td>15.7407</td><td>0.1378</td><td>0.7967</td><td>0.1132</td></tr><tr><td>8.9062</td><td>22.8631</td><td>0.0527</td><td>0.9603</td><td>0.0203</td></tr></tbody></table>

Moreover, letting

$$ \tilde{\psi}(z) = \frac{2}{(1+z)^2} \psi\left(\frac{1-z}{1+z}\right), $$

we get

$$ \int_{0}^{+\infty} \psi(\sigma) d\sigma = \int_{-1}^{1} \tilde{\psi}(z) dz \approx \sum_{i=0}^{N} \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 the Table 1 the weights and the nodes of this quadrature formula for $N = 8$. By identification, we deduce that the quadrature coefficients in (11) 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²⁶ for the approximation of Caputo fractional-order derivative.

# 3 | A FIRST NUMERICAL SCHEME

In this section, we first present the numerical scheme. After that, we present different numerical simulations with both Gauss-Laguerre and Gauss-Jacobi quadrature methods. Finally, we compare these results in order to choose the appropriate quadrature that will be used in the second numerical scheme. In the sequel, we take $t \in [0, T]$, $x \in [0, A]$ and we consider periodic boundary conditions, ie, $u|_{x=0} = u|_{x=A}$. Since we run the numerics with an initial data that provides a wave that moves from the left to the right, we expect our computations to be physically relevant until this wave reaches the right boundary.

## 3.1 | Presentation of the model

Let $\Delta t > 0$ be a fixed time step and $h > 0$ be a fixed space step. We set for all $m = 0, \dots, N_p, t_m = m\Delta t$ and for all $j = 1, \dots, K, x_j = jh$. We use the following notations:

$$ \forall i = 1, \dots, N, \quad u^m(x) \approx u(m\Delta t, x), \quad \psi_i^m(x) \approx \psi(m\Delta t, x, \sigma_i). $$

For the PDE in (10), we suggest a semi-implicit scheme, where the half derivative and the nonlinear term are explicitly approximated. We use the Crank-Nicolson approximation in time. Thus, the semi-discretized scheme associated to (10) is given by
$$ \frac{u^{m+1} - u^m}{\Delta t} + \sqrt{\nu} \sum_{i=1}^{N} w_i \left( \frac{2}{\pi} u^m - \sigma_i^2 \psi_i^m \right) = \alpha \left( \frac{u^{m+1} + u^m}{2} \right)_{xx} - \left( \frac{u^{m+1} + u^m}{2} \right)_x - \beta \left( \frac{u^{m+1} + u^m}{2} \right)_{xxx} - \gamma \left( \frac{(u^m)^2}{2} \right)_x. \quad (13) $$

For the ODE in (10), we use the following implicit-explicit scheme

$$ \frac{\psi_i^{m+1} - \psi_i^m}{\Delta t} = -\sigma_i^2 \psi_i^{m+1} + \frac{2}{\pi} u^m. \quad (14) $$

Then, we apply the Fourier transform in space to (13) and (14). We get

$$ \hat{u}^{m+1} \left(1 + \frac{\Delta t}{2}(\alpha\xi^2 + i\xi - i\xi^3\beta)\right) = \hat{u}^m \left(1 - \frac{\Delta t}{2}(\alpha\xi^2 + i\xi - i\xi^3\beta)\right) \\ - \sqrt{\nu}\Delta t \sum_{i=1}^{N} w_i \left(\frac{2}{\pi}\hat{u}^m - \sigma_i^2 \hat{\psi}_i^m\right) - i\gamma\Delta t \xi \frac{\widehat{(u^m)^2}}{2} \quad (15) $$

and

$$ \hat{\psi}_i^{m+1} (1 + \sigma_i^2 \Delta t) = \hat{\psi}_i^m + \frac{2}{\pi} \Delta t \hat{u}^m, \quad (16) $$

where the initial conditions are approximated as follows:

$$ \forall i = 1, \dots, N, \quad \hat{\psi}_i^0 = 0, \quad \hat{u}^0 \approx \hat{u}_0. $$

Hence, our scheme is a 2-step scheme: Let the approximated solution $(u^m, \psi^m)$ be given, we first calculate an approximate solution $\psi^{m+1}$ of (16). Then, we use this solution to approximate the solution $u^{m+1}$ of (15), and so on.

## 3.2 | Numerical results and discussion

Similar to Chen et al³ and Manoubi,⁵ we choose an initial datum $u_0$, which provides a small amplitude and long-wave Korteweg-de Vries (KdV) soliton for $\alpha = \nu = 0$, $\beta = 1$, and $\gamma = 6$. We denote by $x_0$ the middle of the space interval, then the initial datum reads

$$ u_0(x) = 0.32 * \operatorname{sech}^2(0.4 * (x - x_0)). \quad (17) $$

All numerical simulations are realized with the time step $\Delta t = 0.1$, the space step $h = 0.1$, and on the space interval $[0, A] = [0, 500]$. In a first step, we study the convergence of the numerical schemes (15) and (16) using the quadrature coefficients of Gauss-Laguerre and Gauss-Jacobi. In a second step, we compare these numerical solutions with those constructed using the scheme (68) developed in Manoubi.⁵ In this paragraph, we take the following parameters values $\alpha = \nu = 1$, $\gamma = \beta = 0$, and $T = 100$, which correspond to the nonlinear case without geometric dispersion. In Figure 1, we present the numerical solution associated to the schemes (15) and (16) increasing the number of Gauss-Laguerre quadrature nodes; more precisely, we take $N = 30, N = 40, N = 50, N = 60$, and $N = 80$. Moreover, we present the $L^2$ and $L^\infty$ decay rates of the solution with $N = 40, N = 50, N = 60$, and $N = 80$. We observe that the numerical solutions converge for $N \ge 60$. However, we do not observe the qualitative convergence of the decay rates even for $N = 80$, which represents an important number of points. Also, we observed this slow convergence using other parameters values of the scheme. This may be due to the poor quality of the algorithm of the Gauss-Laguerre quadrature; this fact has been observed by other authors.²⁷,²⁸ Later, this method was mathematically analyzed by Diethelm.²⁶ Diethelm identified that this unsatisfactory behavior is due to the asymptotic characteristics of the integrand function at the end points of the interval $[0, +\infty)$. On the basis of this result, Diethelm proposed an alternative quadrature formula, namely, the Gauss-Jacobi quadrature, that gives rise to a much better overall performance of the scheme and that is more promising for integrands with an algebraic decay.

In Figure 2, we present the approximate solution of the schemes (15) and (16) with the Gauss-Jacobi quadrature coefficients and the $L^2$ and $L^\infty$ decay rates of the solution. We observe that the numerical solutions and the decay rates converge for a number of nodes greater than $N = 15$. We have constructed different numerical simulations using different parameters values, and we get the same observation.
<figure><img src="image_4.png" /><figcaption>**FIGURE 1** Numerical solution and decay rates using Gauss-Laguerre rule for different *N* when $\alpha = \nu = 1$, $\beta = \gamma = 0$ [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

In Figure 3, we compare the numerical solutions and decay rates resulting from the schemes (15) and (16) using the Gauss-Laguerre and Gauss-Jacobi rules with those arising from the numerical scheme (68) developed in Manoubi.⁵ We use 80 points for Gauss-Laguerre quadrature and 15 points for Gauss-Jacobi quadrature. We observe that the 3 numerical solutions overlap each other. Moreover, it is sufficient to use only *N* = 15 points with the Gauss-Jacobi rule to obtain satisfying numerical results for the numerical solution. However, using Gauss-Laguerre quadrature, we need to use *N* = 50 points or more if we need to get more precise results. Moreover, using Gauss-Jacobi rule, we observe that the decay rates of the solution are smoother than both numerical decay rates resulting from the scheme (68) developed in Manoubi⁵ and from Gauss-Laguerre rule. Taking in account these observations, we choose in the sequel the Gauss-Jacobi quadrature to all numerical simulations. Now, we simulate the numerical schemes (15) and (16) in the nonlinear case. To this end, we take the following parameters values: $\alpha = 2$, $\nu = 0.2$, $\gamma = 1$, $\beta = 0.5$, and $T = 100$.

In Figure 4, we present the numerical solution and $L^2$ and $L^\infty$ decay rates of the schemes (15)and (16) using $N = 15$ points in the Gauss-Jacobi rule. We compare these results with those constructed from the numerical scheme (68) developed in Manoubi.⁵ We observe that the numerical solutions overlap each other. Moreover, the decay rates in $L^2$ and $L^\infty$ norms are very close. Then, we can deduce that both schemes converge. Also, in Table 2, we compare the computation time of solutions in linear and nonlinear cases (in seconds) between the numerical schemes (15) and (16) for different
<figure><img src="image_5.png" /><figcaption>**FIGURE 2** Numerical solution and decay rates using the Gauss-Jacobi for different N with $\alpha = \nu = 1$, $\beta = \gamma = 0$ [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

Gauss quadrature rules and the numerical scheme (68) developed in Manoubi.⁵ We observe that the numerical scheme (68) developed in Manoubi⁵ is relatively expensive in computation time. This is due to the numerical method used to approximate the half-order derivative.

## 4 | A SECOND NUMERICAL SCHEME

To improve the precision and the efficiency of the numerical scheme used before, we construct in this subsection a second numerical scheme associated to (8) following the method proposed in Lombard and Mercier.²⁹
<figure><img src="image_6.png" /><figcaption>FIGURE 3 Comparing solutions and decay rates for different quadrature rule with the scheme (68) of Manoubi⁵ when $\alpha = \nu = 1$, $\alpha = \beta = 0$ [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

## 4.1 | Presentation of the model

The PDE-ODE system given by (9) can be written as

$$ \left\{ \begin{array}{l} \partial_t u(t,x) + \partial_x \left( u + \frac{\gamma}{2} u^2 \right) = -\sqrt{\nu} \int_0^\infty \left( \frac{2}{\pi} u(t,x) - \sigma^2 \psi(t,x,\sigma) \right) d\sigma \\ \qquad + \alpha u_{xx}(t,x) - \beta u_{xxx}(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{array} \right. \quad (18) $$

supplemented with the initial conditions

$$ \begin{array}{l} \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{array} $$
<figure><img src="image_7.png" /><figcaption>FIGURE 4 Comparing solutions and decay rates using Gauss-Jacobi quadrature with the scheme (68) of Manoubi⁵ when $\alpha = 2$, $\nu = 0.5$, $\gamma = 0.5$, $\beta = 0.1$ [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

To approximate the integral in (18), we use the Gauss-Jacobi quadrature with $N$ nodes. We denote by $w_i$ and $\sigma_i$ for $i = 1, \dots, N$ the weights and abscissas, respectively. Hence, the system (18) is written as a first-order system as follows:

$$ \left\{ \begin{aligned} \partial_t u(t,x) + \partial_x(u + \frac{\gamma}{2}u^2) &= -\sqrt{\nu} \sum_{i=1}^{N} w_i \left( \frac{2}{\pi} u(t,x) - \sigma_i^2 \psi_i(t,x) \right) \\ &\quad + \alpha u_{xx}(t,x) - \beta u_{xxx}(t,x), \quad t > 0, x \in \mathbb{R}, \\ \partial_t \psi_i(t,x) &= -\sigma_i^2 \psi_i(t,x) + \frac{2}{\pi} u(t,x), \quad t > 0, x \in \mathbb{R}, \quad i = 1, \dots, N, \end{aligned} \right. \quad (19) $$

endowed with the initial conditions

$$ \begin{array}{l} \forall x \in \mathbb{R}, \forall i = 1, \dots, N, \quad \psi_i(0,x) = 0, \\ \forall x \in \mathbb{R}, u(0,x) = u_0(x). \end{array} $$
<table><caption>**TABLE 2** Computation time of the numerical solution</caption><thead><tr><th></th><th>Linear case</th><th>Nonlinear case</th></tr></thead><tbody><tr><td>Scheme (68) in Manoubi<sup>5</sup></td><td>104.26</td><td>125.81</td></tr><tr><td>Gauss-Laguerre (N = 80)</td><td>61.195</td><td>61.45</td></tr><tr><td>Gauss-Jacobi (N = 15)</td><td>40.44</td><td>40.24</td></tr></tbody></table>

Now, we note by

$$U = (u, \psi_1, \dots, \psi_N)^T,$$

the vector of (N + 1) unknowns, by

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

and finally

$$S = \begin{pmatrix} -\sqrt{\nu} \sum_{i=1}^{N} w_i \sqrt{\nu w_1} \sigma_1^2 \dots \sqrt{\nu w_N} \sigma_N^2 \\ \frac{2}{\pi} \quad -\sigma_1^2 \quad 0 \quad \dots \quad 0 \\ \frac{2}{\pi} \quad 0 \quad -\sigma_2^2 \quad \dots \quad 0 \\ \vdots \quad \vdots \quad \vdots \quad \ddots \quad \vdots \\ \frac{2}{\pi} \quad 0 \quad \dots \quad \dots \quad -\sigma_N^2 \end{pmatrix}.$$

Thus, the problem (19) is written in the following form:

$$\partial_t U + \partial_x F(U) = S(U) + G_1 \partial_x^2 U + G_2 \partial_x^3 U, \quad (20)$$

where $G_1$ and $G_2$ are diagonal matrices of order $N + 1$. In the sequel, we introduce the so-called splitting scheme.

## 4.2 | The splitting method

Let $\Delta t > 0$, for all $m = 0, \dots, N_p$, we note by $t^m = m\Delta t$ and

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

From Equation 20, we consider the propagation equation

$$\partial_t U + \partial_x F(U) = G_1 \partial_x^2 U + G_2 \partial_x^3 U, \quad (21)$$

and the diffusive equation

$$\partial_t U = S(U). \quad (22)$$

We note by $H_a$ (respectively, $H_b$) the discrete operator of the solution of (21) (respectively, the solution of Equation 22). Then, a Strang splitting method$^{30,31}$ of order 2 between $t_n$ and $t_{n+1}$ is used to solve respectively (21) and (22) as follows:

$$\begin{align} U^{(1)} &= H_b\left(\frac{\Delta t}{2}\right)U^m, \notag \\ U^{(2)} &= H_a(\Delta t)U^{(1)}, \tag{23} \end{align}$$

$$U^{m+1} = H_b\left(\frac{\Delta}{2}\right)U^{(2)}.$$

Here, the constructed operators $H_a$ and $H_b$ are stable and of order 2. Then, the scheme (23) provides an approximation of order 2 in time to the problem (21). In the sequel, we present the discretization of (21) and (22).

### 4.2.1 | Propagation equation 21

Here, $u$ is a solution of the KdV equation. The discretization in time of this equation is provided by a semi-discrete in time scheme, using the Crank-Nicolson scheme for the linear part and Adams-Bashforth scheme (see Kalisch and Bona$^{32}$) for
the nonlinear part. The space approximation of the solutions are performed by standard Fourier methods. First, we note that the approximate solution $\hat{u}^1$ is provided by a fixed-point method that verifies the semi-discrete scheme (of order 2):

$$ \frac{\hat{u}^1 - \hat{u}^0}{\Delta t} = \frac{\hat{u}^1 + \hat{u}^0}{2} (-\alpha \xi^2 - i\xi + i\beta \xi^3) - \frac{i\gamma \xi (\widehat{u^1})^2 + (\widehat{u^0})^2}{4}. \quad (24) $$

Then, for $m = 1, \dots, N_p$, we consider the following semi-discrete scheme:

$$ \frac{\hat{u}^{m+1} - \hat{u}^m}{\Delta t} = \frac{\hat{u}^{m+1} + \hat{u}^m}{2} (-\alpha \xi^2 - i\xi + i\beta \xi^3) - \frac{i\gamma \xi (3(\widehat{u^m})^2 - (\widehat{u^{m-1}})^2)}{4}. \quad (25) $$

This scheme has local truncation error of order $(\Delta t)^2$, and a second-order convergence is observed (for more details, see Kalisch and Bona³²).

### 4.2.2 | Diffusion equation 22

We can solve mathematically (22) as follows:

$$ \mathbf{H}_b(\frac{\Delta t}{2})U = e^{\mathbf{S} \frac{\Delta t}{2}} U. \quad (26) $$

In addition, for all $N > 0$, the exponential of the matrix $\mathbf{S} \frac{\Delta t}{2}$ is calculated numerically using the “scaling and squaring” method with a (6, 6) Padé approximation³³ that corresponds, with Matlab, to the function “expm($\mathbf{S} \frac{\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 (23) is a 3-step scheme. Since the scheme (25) as well as the Strang splitting method are of order 2 in time, it follows that (23) is of order 2 in time.

## 4.3 | Numerical results and discussion

In this paragraph, we present some numerical results from the numerical scheme (23) with similar parameters and initial datum as in the case of the first schemes 15 and (16). We analyze the convergence in time, and we study the effects of the number of quadrature points on the numerical solution. To this end, we choose $\alpha = 2$, $\nu = 1$, $\gamma = 3$, $\beta = 1$, $T = 100$ and $h = 0.05$.

### 4.3.1 | Convergence in time

To study the convergence in time of the numerical solution of (23), we take $N = 40$ Gauss-Jacobi quadrature points. In Figure 5, we present different numerical solutions for $\Delta t = 0.05$, $\Delta t = 0.1$, $\Delta t = 0.2$, $\Delta t = 0.3$, and $\Delta t = 0.4$, as well as the decay rates in norms $L^2$ and $L^\infty$ in this case. We observe that the difference between different numerical solutions and decay rates are negligible.

To evaluate the error due to the time discretization of the scheme (23), we denote by $u^m_{\text{Ref}}$ the reference solution, which is the numerical solution for $N = 40$ and $\Delta t = 0.025$. We denote by $E^m(\Delta t) = \|u^m_{\text{Ref}} - u^m\|_2$ the $L^2$ error due to the time discretization of the solution in norm $L^2$. We denote by $u^m$ the numerical solution for different time steps $\Delta t$. In Figure 6,
<figure><img src="image_8.png" /><figcaption>**FIGURE 5** Comparison of the splitting scheme solutions and decay rates in $L^\infty$ and $L^2$ norms where $\alpha = 2$, $\nu = 1$, $\gamma = 3$, $\beta = 1$ for different time steps $\Delta t$ [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

we illustrate the $L^2$-error $E^m$ in terms of the time step $\Delta t$ in the log-log scale. We observe that the measured values are close to a straight line with slope 2. This numerically confirms that the numerical scheme (23) is of order 2 in time.

### 4.3.2 | Influence of the number of Gauss-Jacobi quadrature points

To analyze the influence of the number of the grid points $N$ on the numerical solution of (23), we denote by $u_{\text{Ref}}^m$ the reference solution, which is the numerical solution for $\Delta t = 0.05$ and $N = 40$. We denote by $u_N^m$ the numerical solution for different values of $N$. We define the error $E_{Sp}^m(N) = \|u_{\text{Ref}}^m - u_N^m\|_2$ in norm $L^2$. In Figure 7, we present $E_{Sp}^m$ in terms of $N$. We observe that as expected the error decreases as the number of $N$ increases. Also, the measured values are close to a straight line with slope -5, which let us conclude that $E_{Sp}^m \approx C(\frac{1}{N})^5$ where $C$ is a constant.
<figure><img src="image_9.png" /><figcaption>FIGURE 6 Error of the time discretization [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 7 Error of the splitting scheme in terms of N [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

5 | DISPERSION ANALYSIS

We recall the definition of the half-order derivative of Riemann-Liouville

$$D^{1/2}u(t) = \frac{1}{\sqrt{\pi}} \partial_t \int_0^{+\infty} \frac{u(s)}{\sqrt{t-s}} ds. \qquad (27)$$

In this section, we present the dispersion analysis of the model (2) in the linear case ($\gamma = 0$) using 2 methods. First, let us define the Laplace and Fourier transforms in time and space, respectively.

$$F_t(f)(\tau) = \tilde{f}(\tau) = \int_0^{+\infty} f(t)e^{-\tau t} dt,$$

$$F_x(f)(k) = \hat{f}(k) = \int_{-\infty}^{+\infty} f(x)e^{-ikx} dx,$$

where $\tau$ is the angular frequency and $k$ is the wavenumber.

The Laplace-Fourier transform of a function $f$ reads

$$\hat{f}(\tau, k) = \int_{0}^{+\infty} e^{-t\tau} \left( \int_{-\infty}^{+\infty} e^{-ixk} f(t, x) dx \right) dt.$$

The first method consists in applying Laplace-Fourier transform to (2) in order to derive the dispersion relation between $\tau$ and $k$. Then, we define the symbol of the half-order derivative (27)
$$ \chi(\tau) = \sqrt{i\tau} = \sqrt{\tau} \left( \frac{1}{\sqrt{2}} + i \frac{1}{\sqrt{2}} \right). $$

We get from (2) the following identity:

$$ (i\tau + \sqrt{\nu}\sqrt{i\tau})\hat{u} + (ik - i\beta k^3)\hat{u} = -\alpha k^2 \hat{u}. $$

We set

$$ \alpha(\tau, k) = -i\beta k^3 + \alpha k^2 + ik + i\tau. $$

Hence, the dispersion relation is given by

$$ \alpha(\tau, k) + \sqrt{\nu}\chi(\tau) = 0. \qquad (28) $$

We note that we get the same dispersion relation as in Chen et al³ where the authors studied the linear dispersion of (2) when the nonlocal diffusive-dispersive term is described by the Caputo half-order derivative with an initial condition $u(0) = 0$.

The second method is devoted to the linear dispersion relation of the system (20) (we take $\gamma = 0$). To this end, we apply Fourier transforms in time and space to (20), we get a set of ($N + 1$) equations whose determinant must be vanishing. Then, the matrix form of the resulting system is given by

$$ A\hat{U} = 0, $$

where

$$ A = \begin{pmatrix} \alpha(\tau, k) + \frac{2}{\pi}\sqrt{\nu}\sum_{i=1}^{N}w_i - \sqrt{\nu}w_1\sigma_1^2 & \cdots & \cdots -\sqrt{\nu}w_N\sigma_N^2 \\ \frac{2}{\pi} & -\sigma_1^2 - i\tau & 0 & \cdots & 0 \\ \frac{2}{\pi} & 0 & -\sigma_2^2 - i\tau & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{2}{\pi} & 0 & \cdots & \cdots & -\sigma_N^2 - i\tau \end{pmatrix}. $$

Thus, the dispersion relation between $\tau$ and $k$ is given by

$$ \alpha(\tau, k) + \frac{2}{\pi}\sqrt{\nu}\left(\sum_{j=1}^{N}w_j - \sum_{j=1}^{N}\frac{w_j\sigma_j^2}{\sigma_j^2 + i\tau}\right) = 0. \qquad (29) $$

We denote by

$$ \tilde{\chi}(\tau) = \frac{2}{\pi} \left( \sum_{j=1}^{N} w_j - \sum_{j=1}^{N} \frac{w_j \sigma_j^2}{\sigma_j^2 + i\tau} \right), $$

then the dispersion relation (29) can be written in a similar way as (28)

$$ \alpha(\tau, k) + \sqrt{\nu}\tilde{\chi}(\tau) = 0. \qquad (30) $$

It is worth to note that when we consider the main Equation 2 without nonlocal diffusive-dispersive term, ie, when $\nu = 0$, the dispersion relations (28) and (30) reduced to the KdV-dispersion relation, namely, $\alpha(\tau, k) = 0$.

# 6 | CONCLUSION

In this article, we propose different numerical schemes to approximate the solution to an asymptotical water wave model where the nonlocal viscosity is described by the Riemann-Liouville half-order derivative. We compare different numerical results of different numerical schemes where the half derivative is described by a diffusive representation. Moreover, we compare our numerical results with those given in Minoubi.⁵
ACKNOWLEDGEMENTS

This work was supported by PHC Utique ASEO and exchange program CNRS/DGRS “modèles non locaux en dynamique des fluides.” A part of this article was performed during the stay of the second author in the LAMFA with the support of the exchange program “SSHN2015.”

ORCID

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

REFERENCES

1. Liu P, Orfila A. Viscous effects on transient long wave propagation. *J Fluid Mech.* 2004;520:83-92.

2. Dutykh D, Dias F. Viscous potential free surface flows in a fluid layer of finite depth. C.R.A.S, Série I. 2007;347:113-118.

3. Chen M, Dumont S, Dupaigne L, Goubet O. Decay of solutions to a water wave model with a nonlocal viscous dispersive term. *Discrete Contin Dyn Syst.* 2010;27:1473-1492.

4. Dumont S, Duval J-B.. Numerical investigation of asymptotical properties of solutions to models for water waves with non local viscosity. *Int J Numer Anal Model.* 2013;10(2):333-349.

5. Manoubi I. 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 Contin Dyn Syst.* 2014;19:2837-2863.

6. Galucio A, Deü J.-F, Dubois F. The $G^a$-scheme for approximation of fractional derivatives: application to the dynamics of dissipative systems. *J Vib Control.* 2008;14:1597-1605.

7. Galucio A, Deü J, Mengué S, Dubois F. An adaptation of the Gear scheme for fractional derivatives. *Comput Methods Appl Mech Eng.* 2006;195:6073-6085.

8. Koh C, Kelly J. Application of fractional derivatives to seismic analysis of base-isolated models. *Earthq Eng Struct Dyn.* 1990;19:229-241.

9. Li J-R. A fast time stepping method for evaluating fractional integrals. *SIAM J Sci Comput.* 2010;31:4696-4714.

10. Shokooh A, Suarez L.. A comparison of numerical methods applied to a fractional model of damping materials. *J Vib Control.* 1999;5:331-354.

11. Zhang W, Shimizu N. Numerical algorithm for dynamic problems involving fractional operators. *JSME Int J Ser C.* 1998;41:364-370.

12. Zhang J, Xu C. Finite difference/spectral approximations to a water wave model with a nonlocal viscous term. *Appl Math Model.* 2014;38:4912-4925.

13. Liu Y, Du Y, He S, Gao W. Finite difference/finite element method for a nonlinear time-fractional fourth-order reaction-diffusion problem. *Comput Mathe Appl.* 2015;70:573-591.

14. Montseny G. *Représentation Diffusive*. Hermes Science Publications; 2005.

15. Montseny G. Diffusion monodimensionnelle et intégration d'ordre 1/2. *Internal LAAS Report N.* 91232; 1991.

16. Montseny G, Audounet J, Mbodge B. Modèle simple d'amortisseur viscoélastique. Application à une corde vibrante. In: RF Curtain, A Bensoussan, eds. *Lecture notes in Control and Information Sciences*, Vol. 185. Washington: JL.Lions-Springer Verlag; 1993:436-446.

17. Montseny G, Audounet J, Mbodge B. Optimal models of fractional integrators and application to systems with fading memory. *IEEE Int Conf Syst, Man Cybern, Le Touquet Fr.* 1993;17-20:65-70.

18. Staffans O. Well-posedness and stabilizability of a viscoelastic equation in energy space. *Trans Am Math Soc.* 1994;345:527-575.

19. Audounet J, Giovangigli V, Roquejoffre J. A threshold phenomenon in the propagation of a point-source initiated flame. *Phys D.* 1998;121:295-316.

20. Helie T, Matignon D. Diffusive representations for the analysis and simulation of flared acoustic pipes with visco-thermal losse. *Math Models Methods Appl Sci.* 2006;16:503-536.

21. Matignon D, Prieur C. Asymptotic stability of linear conservative systems when coupled with diffusive systems. *ESAIM: Control Optim Calc Var.* 2005;11:487-507.

22. Montseny G, Audounet J, Matignon D. Diffusive representation for pseudo-differentially damped nonlinear systems. *Nonlinear control in the year 2000.* 2000;2:163-182.

23. Abramowitz M, Stegun I. Handbook of mathematical functions with formulas, graphs and mathematical tables. *NBS Applied Mathematics Series 55, National Bureau of Standards*; 1964; Washington, D.C.

24. Yuan L, Agrawal O. A numerical scheme for dynamics systems containing fractional derivatives. *J Vib Acoust.* 2002;124:321-324.

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

26. Diethelm K. An investigation of some nonclassical methods for the numerical approximation of Caputo-type fractional derivatives. *Numer Algorithms.* 2008;47:361-390.

27. Lu J.-F, Hanyga A. Wave field simulation for heterogeneous porous media with singular memory drag force. *J Comput Phys.* 2005;208:651-674.

28. Schmidt A, Gaul L. On a critique of a numerical scheme for the calculation of fractionally damped dynamical systems. *Mech Res Commun.* 2006;33:99-107.
29. Lombard B, Mercier J-F. Numerical modeling of nonlinear acoustic waves in a tube connected with Helmholtz. *J Comput Phys.* 2014;259:421-443.

30. Holden H, Karlsen K, Risebro N, Tao T. Operator splitting for the KDV equation. *Math Comput.* 2011;80:821-846.

31. LeVeque R. *Numerical Methods for Conservation Laws*. 2nd ed. Basel: Birkhäuser-Verlag; 1992.

32. Kalisch H, Bona J. Models for internal waves in deep water. *Discrete Contin Dyn Syst.* 2000;6:1-20.

33. Moler C, Van Loan C. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. *SIAM Rev.* 2003;45:3-49.