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

In this paper, we study the water wave model with a nonlocal viscous term 
 
\begin{equation*} 
 u_t + u_x + \beta u_{x x x} + \frac{\sqrt \nu}{\sqrt \pi}\frac{\partial}{\partial t} \int_0^t\frac{u(s)}{\sqrt{t-s}} ds + u u_x = v u_{xx}, 
\end{equation*} 
 
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. We prove the well-posedness of the equation and we investigate theoretically and numerically the asymptotical behavior of the solutions. Also, we compare our theoretical and numerical results with those given in [4] for a similar equation.


(Communicated by Roger Temam)
Abstract. In this paper, we study the water wave model with a nonlocal viscous term ds is the Riemann-Liouville half derivative. We prove the well-posedness of the equation and we investigate theoretically and numerically the asymptotical behavior of the solutions. Also, we compare our theoretical and numerical results with those given in [4] for a similar equation.

1.
Introduction. The mathematical modeling and analysis of water wave propagation are challenging topics. Intensitive works and progress have been done by a number of authors since L. Euler, J. Boussinesq and D. Korteweg and G. de-Vries for the derivation of the models. 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 [2]. Modeling the effects of viscosity on the propagation of long waves is an important challenge and has been investigated since the time of Stokes and has received a much of interest in the last decade. The exploratory work for this issue is due to T. Kakutani and K. Matsuuchi [14] who have examined the effect of viscosity on the propagation of nonlinear long gravity waves. Besides, P. Liu and T. Orfila [15], and D. Dutykh and F. Dias [8] 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 the linear 3D case. A one-dimensional nonlinear system is presented in [7]. Computing the decay rate for solutions to dispersive-diffusive equations is also a challenging issue. In their pioneering work [1], C. Amick et al. have studied the large-time behavior of solutions to the initial value problem for the Kortewegde-Vries equation and for the regularized long-wave equation, with a dissipative term appended. More theoretical and numerical methods and details can be found in [3,5,6,12,13] for computing the decay rates for solutions to dissipative PDEs.
In their recent work [4], M. Chen et al. investigated theoretically and numerically the decay rate for solutions to the following water wave model We denote as β, ν and α the parameters dedicated to balance the effects of viscosity and dispersion against nonlinear effects. Particularly, the authors in [4] consider (1) with β = 0 supplemented with the initial condition u 0 ∈ L 1 (R) ∩ L 2 (R). They proved that if u 0 L 1 (R) is small enough, then there exists a unique global solution u ∈ C(R + ; L 2 x (R)) ∩ C 1 (R + ; H −2 x (R)). In addition, u satisfies In this paper, we are interested in the model (1) where the fractional term is described by the Riemann-Liouville half derivative instead of that of Caputo, namely We prove the local and global existence of solutions to problem (3) when β = 0 using a fixed point theorem. Also we study theoretically the decay rate of the solutions in this case. Moreover we perform numerical simulations on the decay rate of the solutions for different values of the parameters α, ν and β. In addition, we are interested to compare our results with those given in [4]. The outline of this article is as follows: in Section 2 we give the main results and notations. In Section 3, we define the solution of an ordinary fractional equation with the Riemann-Liouville operator. Moreover we prove Theorem 2.1 which provides solutions to the linear problem when β = 0. Also, we investigate the asymptotical behavior of the kernels. Then in Section 4 we prove Theorem 2.2 which provides the local and global existence of solutions of the nonlinear problem when β = 0. Besides we study the decay estimates of these solutions. Finally in Section 5, we present the numerical schemes and we perform several numerical simulations for (3). Also, we investigate numerically the decay rate of the solutions.

2.
Main results and notations. In the following, we outline the main results and we introduce the notations used in the sequel.
2.1. Statement of the theoretical results. In all the following, represents the usual convolution product and is the time-space convolution product defined by v w(t, x) = = max(−x, 0) presents the negative part of x ∈ R. We introduce also the kernels and We begin with the following initial value problem for the linear equation.
where the initial data u 0 ∈ L 2 (R). We have the following result.
Theorem 2.1. For all u 0 ∈ L 2 (R), there exists a unique solution u for the linear problem (6) where K RL is the convolution kernel given by (4). Besides, we have ). Now we state a global well-posedness result for the following nonlinear equation supplemented with a small initial data u 0 .
Theorem 2.2. Let u 0 ∈ L 2 (R), then there exists a unique local solution u ∈ C([0, T ); L 2 x (R)) of (8). Moreover for u 0 ∈ L 1 (R) ∩ L 2 (R), there exists a positive constant C 0 > 0 that depends on u 0 such that if u 0 L 1 (R) is small enough, there exists a unique global solution u ∈ C(R + ; L 2 x (R)) ∩ C 1/2 (R + ; H −2 x (R)) of (8) given by where K RL and N are given by (4) and (5), respectively. In addition, we have the following estimate Remark 1. In other words (10) provides the decay estimate for solutions.

Notations and definitions.
Let us introduce some notations that we will use in the sequel. A function f is of exponential order α if there exists a constant M > 0 and t 0 ≥ 0 such that . Now let f be a piecewise continuous causal function of exponential order α. We define the Laplace transform of f by for all τ such that Re(τ ) > α.

IMEN MANOUBI
The Fourier transform of a function f ∈ L 1 (R) is defined as follows We denote by L p (R) the usual Lebesgue space and by H s (R) with s ∈ R the usual Sobolev space. Also we stand by C k+α (R) the usual Hölder space with k ∈ N and α ∈ (0, 1).
3. Proof of Theorem 2.1. In this section, we introduce a fractional differential equation (FDE) where the fractional term is described by the Riemann-Liouville half derivative and we prove the well-posedness of the associated initial value problem. Then we investigate the linear homogeneous problem given by (3.2). Finally, we state results corresponding to the linear inhomogeneous problem. To begin with, we recall some preliminary results.

3.1.
Auxiliary results. Let f : R + → R be a locally integrable function and consider the following integro-differential equation supplemented with the initial condition v(0) = v 0 ∈ R. We seek a solution for this equation. To this end , we integrate (11) between 0 and t, so v satisfies the following integral equation Conversely, let v be a solution of (12).  Proposition 1. Let f be a function in C(R + ) and let v 0 ∈ R. The problem (11) supplemented with the initial condition v(0) = v 0 admits a unique solution v ∈ C 1/2 (R + ). Moreover, v is given by where Proof of Proposition 1. We begin with proving the existence of a solution to problem (13). To this end, we apply the Laplace transform to equation (11). Since This is equivalent toṽ To continue the proof, we need the following lemma Proof of Lemma 3.2. Applying the Laplace transform to (14), we get Using Fubini-Tonelli theorem, it follows Splitting the integrals and using a change of variables it follows Now using the properties of the Laplace transform, the identity (15) is equivalent to

IMEN MANOUBI
By identification, we conclude that the solution v of (11) is given by (13). Let u and u be two solutions of (13) in C(R + ) supplemented with the same initial condition u 0 . Let w = u −ũ so w(0) = 0 and w satisfies the following equation To complete this part of the proof, we state the following proposition.
Proposition 2. Let λ ∈ C and w a continuous function on R + which is solution of the following integral equation where w(0) = 0. Hence, we have the following result Proof of Proposition 2. Let w be a solution of (17), so we have the following inequality We denote by It follows that Let T > 0 be small enough such that 2 √ π √ T + |λ|T < 1, then m(t) = 0 for all t ∈ [0, T ]. This result may be extended for all T ≥ 0 by proving recursively m(kT ) = 0 for all k ∈ N and k ≥ 2.
Hence the Proposition 2 is proved.
Using Proposition 2 we conclude that the solution of equation (16) is identically 0 on R + and thus the problem (12) admits a unique solution in C(R + ).
In the following, we study the regularity of the solution v of the equation (12). To this end, we rewrite the function v in (13) as By an integration by parts, we see that

Remark 2.
If v 0 = 0 then the solution v of the equation (11) belongs to C 1 (R + ).

3.2.
The linear homogeneous problem. In this subsection, our aim is to prove Theorem 2.1. To this end, we follow the method introduced by M. Chen et al. in [4]. This method consists on applying the so-called Laplace-Fourier transform to the integro-differential equation. Then we use the inverse Fourier transform and the inverse Laplace transform to get the solution in (t, x) variables. We start by providing some estimates on the convolution kernel.

3.2.1.
Estimates on the convolution kernel. First we recall the expression of the convolution kernel Proposition 3. The kernel K RL is non negative Proof. In fact we observe that This implies that the kernel K RL is non negative.
We have the following estimates on the kernel in L ∞ , L 1 and L 2 norms.
Proof. First, we present a general estimate on K RL for all positive time t. Then we give special estimates for large time t > 1. Hence, for all t ≥ 0 we have We conclude that and that Using Hölder inequality, we get These estimates are given for all time t ≥ 0 and are valid if 0 ≤ t ≤ 1. We note that these results are similar to those given in [4] for the estimates of the kernel with the Caputo half derivative. Now, consider t > 1 and rewriting K RL as follows

IMEN MANOUBI
This implies that Thus we conclude that Using Hölder inequality, we get This completes the proof.

3.2.2.
Computing the solution. We formally apply the Laplace-Fourier transform to (6), we get We denote byû the Laplace-Fourier transform of u. For the state of simplicity, we also denote byû 0 the usual Fourier transform of u 0 . Hence, the above equation is equivalent to Lemma 3.4. The kernelK RL defined by (20) is the Laplace-Fourier transform of Proof. We note thatK RL is the Fourier transform in space of It follows that In the sequel, we need the following proposition Hence, we begin with providing the inverse Laplace transform in time of the following functionf where χ is the characteristic function. Multiplying (21) by e |x|/2 e −x − , and using the properties of the Laplace transform it follows Thus, we deduce from Proposition 4 that Hence An integration by parts of the right hand side of (22) followed by a change of variables yields Finally we deduce that the kernel is given by This completes the proof of the Lemma.
Recalling (19) and using the inverse Fourier transform we get Hence by identification we conclude that This result is given for an initial data u 0 ∈ L 1 (R) ∩ L 2 (R). We can extend it straightforwardly for all u 0 ∈ L 2 (R) by a density argument.

3.2.4.
Regularity of the solution. In this paragraph we first prove that the solution u ∈ C(R + ; L 2 x (R)). Then we show that u ∈ C 1/2 (R + ; H −2 x (R)). Finally we check that u ∈ C 1 (]0, +∞[; H −2 x (R)). Let t 0 > 0 and t such that t → t 0 . Using the properties of the convolution product in L p -spaces we get . On the other hand the function t → K RL (t, x) is continuous for all t > 0. Then using the Lebesgue dominated convergence theorem we conclude that for all t > 0 we have lim To prove the continuity on 0 of the solution, we rewrite the kernel K RL as follows The proof is based on the fact that K 1 (t, ·) is an approximate of the identity in L 2 x (R), which implies that lim Moreover using the Lebesgue dominated convergence theorem, we deduce that which completes the first part of the proof.

3.3.
The linear inhomogeneous problem. This subsection is devoted to the following problem supplemented with the initial condition u(0, x) = 0 for all x ∈ R. It is worth to note that since u(0, x) = 0, the Riemann-Liouville half derivative is equal to the Caputo half derivative, namely Thus we recall the results given in [4], Subsection 2.3.

Proposition 5.
Consider a function f ∈ C(R + ; L 1 x (R)), then there exists a unique solution u to (28) given by where N is given by (5). Moreover we have Lemma 3.5. There exists a constant C > 0 such that for all t > 0, we have the following estimates on the kernel

Proof of Theorem 2.2. Consider the nonlinear fractional PDE
supplemented with the initial condition We can obtain, at least formally, that a solution of (29)-(30) is a solution of the following (mild) Duhamel formulation We solve (31) by applying a contraction principle using a Picard fixed point theorem to a functional Φ defined by in an appropriate metric space that we will define in the sequel.

Local existence.
In this paragraph, we prove the existence of a local solution to (29) with an initial condition in L 2 (R). Let T > 0 and set continuous } X T is a Banach space when endowed with the norm u X T := sup t∈[0,T ] u(t) L 2 x (R) . In the following we give some properties satisfied by the functional Φ on X T with u 0 ∈ L 2 (R). Let u ∈ X T ,we have

This implies that
Moreover we have Using Lemma 3.5, we obtain Using the fact that for all t ≥ 0 we have By (33) and (34) we deduce that for all u ∈ X T Let u and v in X T then for all t > 0 we have Using Hölder inequality we get This implies that Also using Theorem 2.1 and Proposition 5, we show that if u 0 ∈ L 2 (R) then Φ is well defined, i.e. for all t ≥ 0 the map t → Φ(u)(t) is continuous for all u ∈ X T . Next we define a set B 0 invariant under the action of Φ. Let u 0 ∈ L 2 (R) and take R 0 = 2C 1 u 0 L 2 . Let B 0 be the closed ball in X T of radius R 0 centered at the origin. Thanks to (35) and (36) we get and In conclusion, we take T > 0 small such that C R0 = max(C 2 , C 3 )R 0 T 1/4 ≤ 1 2 . Hence with this choice we get Φ(B 0 ) ⊂ B 0 , and thus the map Φ is a contraction on B 0 . Using the fixed point theorem we deduce that there exists a unique fixed point u of the functional Φ in the ball B 0 .
In the sequel we prove the global existence of the solution.

Global existence.
In the following, we assume that u 0 ∈ L 1 (R) L 2 (R). We define the Banach space endowed with the norm v X = sup t>0 t 1/4 v(t, ·) L 2 x (R) . Now we give some properties satisfied by Φ on X when u 0 ∈ L 1 (R) L 2 (R). First we have for all t ≥ 0 In addition, let u ∈ X then for all t ≥ 0 we have Taking into account that Hence, recalling (39) and (40), we get for all u ∈ X Let u and v be in B and t ≥ 0, we have In addition using Theorem 2.1 and Proposition 5 we show that if u 0 ∈ L 1 (R) L 2 (R) and v ∈ X then Φ(u) ∈ X. Next let u 0 ∈ L 1 (R) L 2 (R) and take R = 2C 1 u 0 L 1 (R) .
We denote by B = {u ∈ X; u X ≤ R} the closed ball of X of radius R centered at the origin. Recalling (41) and (42) we get and In conclusion, we take u 0 L 1 (R) small enough such that C R = max(C 2 , C 3 )R ≤ 1 2 . Hence with this choice we get and thus the map Φ is a contraction on B. We deduce that Φ admits a unique fixed point in B given by (31).
For the regularity of the solution, we note that u is a superposition of the solution of the linear homogeneous problem investigated in Subsection 3.2 and the solution of inhomogeneous problem investigated in Subsection 3.3; this implies that u ∈ C(R + ; L 2 x (R)) ∩ C 1/2 (R + ; H −2 x (R)).

Decay estimates.
In this subsection, we prove the decay rate of solutions for small initial data. To this end, we proceed in two steps.
First step. Let ε > 0 be sufficiently small and let u 0 ∈ L 1 (R) L 2 (R) such that We first consider the case 0 < t < 3, we observe that Hence we conclude that for all 0 < t < 3 We now consider the case t ≥ 2, we have It follows that For all t and s ∈ R + , we define Hence, we get for all t ≥ 2 Using a change of variables on the integrals in the right hand side of (48), we obtain Thus there exists a constant C > 0 such that In addition, Besides using the fact that t ≥ 2 we obtain Collecting these results we conclude that Recalling (50) and (51) we have . We define f (X) = C X 2 − X + Cε. Hence M 2 satisfies the inequality f (X) ≥ 0. We note that there exists two distinct positive roots of f , namely, Hence given C and C we can choose ε > 0 small enough such that f (2Cε) = Cε(4C Cε − 1) < 0 which implies that M 2 (t) ≤ 2Cε. Hence for t ≥ 2, we have This completes the proof of the upper bound in L 2 norm.
Second step. Let u 0 L 1 (R) ≤ ε with ε > 0 small. Hence for all t ∈ R + we have Let 0 < t < 3, we observe that Hence we deduce that for all 0 < t < 3 It follows that We define for all t ∈ R + M ∞ (t) = sup On the one hand we have Using (56) and the change of variables σ = s/t, it follows On the other hand This implies that Using the change of variables σ = s/t yields

Hence we get
Recalling (55), (57) and (58), we deduce that for all t ≥ 2 Recalling (52), we conclude that CM 2 (t) 2 ≤ Cε and CM 2 (t)M ∞ (t) ≤ 1 2 M ∞ (t). Hence we get for all t ≥ 2 5. Numerical computation for the fractional PDE. This section deals with the numerical solution of the nonlinear equation supplemented with the initial condition u 0 . Here ν and α are non negative parameters. In the following we use two methods to approximate the solution of (61). The first one is detailed in [4] and is based on the kernel representation of the solution.
The second method is based on the approximation of the Riemann-Liouville half derivative using the Gear scheme [6,11].
5.1. The numerical schemes. We reformulate the problem (61) as follows In the sequel, we take x ∈ [0, A] and we work with periodic boundary conditions i.e u |x=0 = u |x=A . The space approximation of the solutions is performed by standard Fourier methods. Since we run the numerics with an initial data that provides us with 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. Let us develop a first numerical scheme associated to equation (62). First we recall that the solution of (62) with initial condition u 0 is given by where N 0 is given by (14). Let ∆t > 0 be fixed time step and set t n = n ∆t, We denote by t k+ 1 2 = t k+1 +t k 2 = (k + 1 2 )∆t. Also let f lin = αu xx − u x − βu xxx the linear part of f and f nl = −γuu x the nonlinear part of f . We define This is equivalent to Subsequently we approximate the integrals in the right hand side of (65) using the following quadrature formula Straightforwardly, this approximation is of order 2 in time. In addition, we have Hence it follows that This approximation is of order 1 in time. We denote by u n the approximate value of u(t n ). Recalling (63), (66) and (67) and then applying the Fourier transform in space to the solution yieldŝ Hence the numerical scheme (68) is of order 1 in time due to the explicit approximation of the nonlinear term. Moreover, this scheme is semi-implicit. We note that, for a fixed initial condition u 0 , we get a unique numerical solution which implies that the scheme (68) is well posed.
We propose now a second scheme based on the Gear operator. For that purpose, we follow the approach proposed in [6]. This method consists on approximating the Riemann-Liouville half derivative using the Gear operator. We give in the following an outline of the Gear scheme developed by A.-C Galucio et al. in [11]. Let G be the Gear operator defined by that approximates the first derivative of u with respect to time. Here δ − is a backward operator defined by (δ − u) n = u n−1 . We formally take the α-derivative of (69) Using Newton binomial formula we get  where (−1) j C j α is given in terms of the Gamma function: . Then the α-derivative of u at each time t n can be approximated by where g j+1 are rational numbers. For illustrative purposes, we present the first five G α -coefficients in Table 1 for three values of α : 1 3 , 1 2 and 3 4 . As it was mentioned in [10], the Gear scheme is a three-level step algorithm, backward in time and second order accurate, for the approximation of classical time derivatives. Also, 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 [11].

5.2.
Numerical results and discussion. Similarly to [4] and [6] we choose an initial datum u 0 which provides a small amplitude and long wave KdV soliton for α = ν = 0, β = 1 and γ = 6. We denote by x 0 the middle of the space interval then Initial data α=0, ν=1 α=5, ν=0 α=5, ν=1 In order to compare the approximated solutions of both schemes (68) and (73) we take an interval of length A = 2500, α = ν = 0.1, β = 0 and γ = 1. We denote by u 1 the numerical solution associated to the scheme (68) and by u 2 the numerical solution associated to the scheme (73). First we take a time step ∆t = 0.1 so we have u 1 L 2 = 0.044 and u 2 L 2 = 0.038 at time T = 50. The maximum difference in L 2 norm between both solutions is equal to u 1 − u 2 L 2 = 0.0066. If we take ∆t = 0.05 we get u 1 L 2 = 0.044 and u 2 L 2 = 0.04. We see that the L 2 norm of solutions are very close. In addition, the difference between both solutions in L 2 norm decreases to u 1 − u 2 L 2 = 0.0047. We deduce that the numerical schemes (68) and (73) are convergent. We note that all numerical simulations presented below are solutions of the scheme (68) with the space step discretization h = 0.2 and the time step discretization ∆t = 0.05. In the following, our purpose is to study the influence of the parameters in (68) on the asymptotical behavior of the numerical solution. Then, we aim to compare these results with those provided in [4].
In Figure 1, we observe the effects of the local and nonlocal viscous terms in the linear case (i.e when γ = 0) for different values of parameters (α, ν) = (0, 1), (α, ν) = (5, 0) and (α, ν) = (5, 1). These values describe respectively the case where we have only nonlocal viscous term, only local viscous term and finally with both terms. The numerical solutions here are provided for T = 100. We note that    we have the same remarks as in the case with the Caputo half derivative. When (α, ν) = (0, 0) the solution of the linear equation presents a travelling wave with speed equal to 1. Comparing this with the case (α, ν) = (0, 1), we see that the nonlocal viscous term slows the wave down significantly and at the same time it enlarges the wave length. In addition, comparing the results when (α, ν) = (0, 0) and (α, ν) = (5, 0), we see that the nonlocal viscous term enlarges the wave length but it does not effect the speed of the wave. Finally, we point out that when both terms are introduced, the profile of the solution is close to that when only nonlocal viscous term is involved. We conclude that the effect of the nonlocal viscosity is more important than that of local viscosity.
In the following we aim to check numerically the theoretical results proved in Theorem 2.1. Also we perform numerical simulations for cases which have not been investigated theoretically. Moreover we compare our results to that presented in [4] with Caputo half derivative. We investigate the decay rate of the solutions of equation (61) in L 2 norm and in L ∞ norm in a large time interval (0, 1000) with different viscous coefficients (ν, α) = (0, 0.1), (ν, α) = (0.1, 1), (ν, α) = (0.1, 0.1), where β = 0 and γ = 0 (linear case) and γ = 1 (nonlinear case). We expect the decay rate of the solution to be as u(t, ·) L 2 x ≈ Ct a or u(t, ·) L ∞ x ≈ Ct a for t large with a, a < 0. Thus we have the following estimates on the ratios It is worth to note that when the nonlocal term is not involved in the simulation, i.e (ν = 0), our numerical results coincide with the one given in [4]. In Figure  2, we simulate the ratio R ∞ for different values of ν and α. One observes that the nonlocal dissipative term produces a larger decay rate compared with the local dissipative term. Also we see that when ν = 0.1, the decay rate of the solution in L ∞ norm is close to (+1). However when ν = 0, the L ∞ -decay rate tends to (+1/2) and this result is expected. These results match well the theoretical results given in Theorem 2.1. In addition, we observe that when ν = 0.1 with α equal to both 0 and 0.1 the curves of decay rates overlap each other. We point out that we have not investigate theoretically the case with only nonlocal term, namely when (ν, α) = (0.1, 0). In Figure 3, we perform the decay rate of the solutions in L 2 norm and the results are similar. We see that when ν = 0 the decay rate approaches to (+3/4) and the curves overlap each other when α is varying. Hence, the influence of α seems to be of smaller than that of ν. However, when ν = 0 the decay rate approaches to (+1/4). Similar numerical simulations of the decay rate of the resulting approximated solutions are performed in Figures 4 and 5 for the nonlinear case (γ = 1) and where the geometric dispersion term is important (β = 6). We point out that there are not theoretical results available in this case. However the results are similar. In Figure 6 we compare numerically the solutions with the Riemann-Liouville half derivative and with the Caputo half derivative at time T = 100. We see that the amplitude of the Riemann-Liouville solution reduces significantly comparing with the Caputo case. However the velocity of the wave is still preserved. This result is expected. In fact, in Theorem 2.2 the decay estimate for the Riemann-Liouville