Splitting methods for the nonlocal Fowler equation

We consider a nonlocal scalar conservation law proposed by Andrew C. Fowler to describe the dynamics of dunes, and we develop a numerical procedure based on splitting methods to approximate its solutions. We begin by proving the convergence of the well-known Lie formula, which is an approximation of the exact solution of order one in time. We next use the split-step Fourier method to approximate the continuous problem using the fast Fourier transform and the finite difference method. Our numerical experiments confirm the theoretical results.


INTRODUCTION
We consider the Fowler equation [11,12]: where u = u(t, x) represents the dune height and I is a nonlocal operator defined as follows: for any Schwartz function ϕ ∈ S(R) and any x ∈ R, We refer to [1,2,6] for theoretical results on this equation.
where F denotes the Fourier transform normalized in (1.11). Thus, I can been seen as a fractional power of order 2/3 of the Laplacian, with the "bad" sign. It will be clear from the analysis below that our results can easily be extended to the case where I is replaced with a Fourier multiplier homogeneous of degree λ ∈]0, 2[, as in [4], and not only λ = 4/3.
Recently, to solve the Fowler equation some numerical experiments have been performed using mainly finite difference approximation schemes [3,4]. However, these schemes are not effective because if we opt for an explicit scheme, numerical stability requires that the time step ∆t is limited by O(∆x 2 ). And, if we choose an implicit scheme, we have to solve a large system which is a computationally expensive operation. Thus, the splitting method becomes an interesting alternative to solve the Fowler model. To our knowledge, there is no convergence result in the literature for the splitting method associated to the Fowler equation. This method is more commonly used to split different physical terms, such as reaction and diffusion terms, see for instance [18]. Splitting methods have also been employed for solving a wide range of nonlinear wave equations. The basic idea of this method is to decompose the original problem into sub-problems and then to approximate the solution of the original problem by solving successively the sub-problems. Various versions of this method have been developed for the nonlinear Schrödinger, Kortewegde-Vries and modified Korteweg-de-Vries equations, see for instance [13,17,19,21]. For the Fowler model (1.1), we consider, separately, the linear Cauchy problem (1.6) ∂v ∂t , and the nonlinear Cauchy problem where ε, η are fixed positive parameters such that ε + η = 1. Equation (1.7) is simply the viscous Burgers' equation. We denote by X t and Y t , respectively, the evolution operator associated with (1.6) and (1.7): where D(t, ·) = F −1 e −t φ I with φ I (ξ) = 4π 2 ηξ 2 − a I |ξ| 4/3 + b I ξ|ξ| 1/3 , and where G is the heat kernel defined by Furthermore, the following L 2 -estimate holds Let us explain the choice of this decomposition. First, we can remark that if we do not consider the nonlinear term in (1.1), the analytical solutions are available using the Fourier transform. Thus, the linear part may be computed efficiently using a fast Fourier transform (hereafter FFT) algorithms. Note also that the Laplacian and the fractional term I cannot be treated separately. Indeed, the equation u t + I[u] = 0 is ill-posed. We next decide to handle the nonlinear term by adding a bit of viscosity in order to avoid shock problems in the standard Burgers' equation. Therefore, the splitting approach presented in this article differs from e.g. the one analyzed in [10], which corresponds to assuming ε = 0 in the above definitions. The splitting operators associated to this approach when I = 0 (which amount to considering alternatively the heat equation and the Burgers equation) have been studied in [15] (as well as other equations involving the Burgers nonlinearity, such as the KdV equation, see also [14]). The main difference with the result presented in [15] is that the operator I is not a differential operator, so its action on nonlinear terms is rather involved , while this point would be needed to compute the Lie commutator between the two operators x v, as in e.g. [15,17]. Note that since we consider the Lie splitting operator, we could in principle use the exact formula established in [7] for the local error. But then again, we face the problem to compute [A, B], which is the only term that we cannot estimate in [7,Theorem 1]. We finally point out that we use smoothing effects associated to the viscous Burgers equation; see Corollaries 3.8 and 3.9.
We motivate this choice by the presence of artificial diffusion in classical numerical schemes used to solve the convection equations. An alternative to reduce this effect is to consider numerical schemes of high order which are usually computationally expensive and do not seem to be very useful for the Fowler model because of the diffusion term.
We consider the Lie formula defined by The alternative definition Z t L = Y t X t could be studied as well, leading to a similar result. Also, the following evolution operators , corresponding to the Strang method [20] could be considered. Following the computations detailed in the present paper for the case (1.10), it would be possible to show that the other Lie formula generates a scheme of order one, and to prove that the Strang method is of order two (for smooth initial data), in the same fashion as in, e.g., [5,17]. This fact is simply illustrated numerically in Section 6, to avoid a lengthy presentation. With Z t L given by (1.10), our main result is: For all u 0 ∈ H 3 (R) and for all T > 0, there exist positive constants c 1 , c 2 and ∆t 0 such that for all ∆t ∈]0, ∆t 0 ] and for all n ∈ N such that 0 n∆t T , Here, c 1 , c 2 and ∆t 0 depend only on T , ρ = max t∈[0,T ] S t u 0 H 2 (R) , and u 0 H 3 (R) . Remark 1.3. It will follow from Lemma 3.11 that for some nonlinear (increasing) function C T depending on T .
In this paper, we begin by estimating the L 2 -stability for error propagation. We next prove that the local error of the Lie formula is an approximation of order two in time. Finally we prove that this evolution operator represents a good approximation, of order one in time, of the evolution operator S t in the sense of Theorem 1.2. Furthermore, we apply Lie and Strang approximations in order to make some numerical simulations using the split-step Fourier experiments.
This paper is organised as follows. In the next section we give some properties related to the kernels G and K, and we prove two fractional Gronwall Lemmas. In Section 3, we establish some estimates on X t , Y t , Z t L and S t . In Section 4, we prove a local L 2 error estimate. Theorem 1.2 is proved in Section 5. We finally perform some numerical experiments which show that the Lie and Strang methods have a convergence rate in O (∆t) and O ∆t 2 ), respectively. Notations.
-We denote by F the Fourier transform of f which is defined by: for all ξ ∈ R, We denote by F −1 its inverse.
-We denote by C T (c 1 , c 2 , · · · ) a generic constant, strictly positive, which depends on parameters c 1 , c 2 , · · · , and T . C is assumed to be a monotone increasing function of its arguments.

PRELIMINARIES
2.1. Properties of the kernels. We begin by recalling the properties of kernels involved in the present analysis.
Proposition 2.1 (Main properties of K, [1]). The kernel K satisfies: For any u 0 ∈ L 2 (R) and t > 0, Proposition 2.2 (Main properties of G, [9]). The kernel G satisfies: x has similar properties to the kernel K. Moreover, for all t > 0, we have Then there exists C T (θ) such that Proof. The proof of this Lemma is well-known and is based on an iteration argument; see for instance [8,Lemma 4.3] or [22,Corollary 2]. We sketch the argument for the sake of completeness. Iterating inequality (2.2) once, we have where β is the beta function. Therefore, we have Iterating the estimate (2.3) n times, with nθ 1, we get the following estimate: with α 0, and whereL T (θ) is a positive constant which depends on T and θ. The lemma then follows from the classical Gronwall Lemma.

Lemma 2.5 (Modified fractional Gronwall Lemma).
Let φ : [0, T ] → R + be a bounded measurable function and P be a polynomial with positive coefficients and no constant term. We assume there exists two positive constants C and θ ∈]0, 1[ such that for all t ∈ [0, T ], Then there exists C T (θ) such that for all t ∈ [0, T ], Proof. Arguing as in the proof of Lemma 2.4, we iterate the previous inequality. After one iteration, we get where we have used the assumptions on P , and Fubini's Theorem again for the last term. Iterating sufficiently many times, we infer like in the proof of Lemma 2.4: Then Using (2.5) to control the second term, and choosing C 1 sufficiently large, we infer: Since This completes the proof.

ESTIMATES ON THE VARIOUS FLOWS
3.1. Estimates on linear flows. In this paragraph, we collect several estimates concerning the convolutions with D, K and G, which will be useful in the estimates of the local error of the scheme.
Proof. For all s ∈ R and all ϕ ∈ H s (R), we have, using (1.3) hence the result.
Then, for all v ∈ H n (R) and all t > 0, There exists C such that for all v ∈ H 2 (R) and all t > 0, Proof. Using Plancherel formula, we have then, from again Plancherel formula, we have hence the first point of the lemma. .
But from the definition of X t , . X s is given by . Thus, using Proposition 3.1 and the first point of this lemma, we get hence the result.
Recalling that K corresponds to D in the case η = 1, we readily infer: For all w ∈ H 2 (R) and all t > 0, where C is a positive constant independent of t and w.
We conclude this paragraph with an analogous result on the heat kernel G: Lemma 3.4. For all w ∈ H 2 (R) and all t > 0, Proof. Proceeding as above, we have: Taking the norm L 2 and using Proposition 2.2, Young's inequality yields hence the result.

3.2.
Estimates on Y t . We now turn to the viscous Burgers' equation (1.7): turns the viscous Burgers' equation into the heat equation: We infer the explicit formula: However, this formula does not seem very helpful in order to establish Proposition 3.6.
In addition, there exists C = C(ε, w 0 H 1 (R) ) such that for all t 0, If in addition w 0 ∈ H n (R), for some n 2, then w ∈ C(R + ; H n (R)) and for all T > 0, Proof. The existence and uniqueness part being standard, we focus on the estimates. The L 2 estimate yields (formally, multiply (3.4) by w and integrate) and the H 1 estimate (differentiate (3.4) with respect to x, multiply by ∂ x w and integrate), The L 2 estimate (3.5) shows that the map t → w(t) 2 L 2 is non-increasing: An integration by parts and Cauchy-Schwarz inequality then yield In order to take advantage of the smoothing effect provided by the viscous part, integrate (3.6) in time and write Gagliardo-Nirenberg inequality yields L 2 , so using (3.7), we infer: In view of Hölder inequality in the last integral in time, where we have used Young inequality ab a 8/7 + b 8 . We infer where we have used (3.7), Hölder inequality and (3.8), successively.
Integrate the H 1 estimate (3.6) with respect to time, and now discard the viscous part whose contribution is non-negative: The first part of the proposition then follows from the Gronwall lemma.
To complete the proof of the proposition, we use the general H n estimate, for n ∈ N: Integrating by parts, the last term is non-positive, since Integrating by parts the first term yields In view of Kato-Ponce estimate [16] (3.10) we have (with f = w and g = ∂ x w) Leaving out the viscous term, Gronwall lemma yields the a priori estimate where C depends only on n ∈ N. In particular, for n = 2, Gronwall lemma implies ∂ 2 x w(t) L 2 w 0 H 2 e C(t 5/8 +t) , where C = C(ε, w 0 H 1 ). We bootstrap, thanks to Gagliardo-Nirenberg inequality again: The last estimates of the proposition then follow from (3.11).
Proof. Differentiating the Duhamel formula (1.8) in space, we have Using Young inequality and inequality (1.9), we infer: In view of Proposition 2.2, this implies: Writing and invoking the fractional Gronwall Lemma 2.4 with θ = 1/4, the lemma follows.
Corollary 3.8. Let n 1 and w ∈ H n (R). Let T, α > 0. If then there exists c depending only on T and α such that Proof. Denote w(t) = Y t w. From Lemma 3.7, where C depends only on α and T . From the proof of the first part of Proposition 3.6 (see the estimate after (3.8)), we infer The corollary then stems from (3.11).

3.3.
Estimates on the splitting operator Z t L . Combining the estimates on X t and Y t established in the previous two sections, we infer: (1) For all u ∈ L 2 (R) and all t > 0, Proof. Multiplying (1.1) by u and integrating with respect to the space variable, we get: because the nonlinear term is zero. Using (1.3) and the fact that u and R (I[u]−∂ 2 xx u)u dx are real, we get We infer where α 0 = − min Re(ψ I ) > 0. The result then follows from the Gronwall lemma.
Lemma 3.11. Let n ∈ N * , u 0 ∈ H n (R) and T > 0. There exists C T ( u 0 H n−1 (R) ) such that the unique mild solution u ∈ C([0, T ]; H n (R)) satisfies Proof. The proof is similar to the one given in Lemma 3.7. Differentiating the Duhamel formula (1.4) in space, we have Using Young inequality and Proposition 2.1, we infer, for any integer n 1: In view of Proposition 2.1, this implies: For n = 1, we use Lemma 3.10 to have The fractional Gronwall Lemma 2.4 with θ = 1/4 then yields where C depends only on T and u 0 L 2 . From (1.9), this implies the lemma in the case n = 1. For n 1, Leibniz rule and Cauchy-Schwarz inequality yield The lemma then easily follows by induction on n.
We will also need the fact that the flow map S t is uniformly Lipschitzean on balls of H 2 (R). Proposition 3.12. Let T, R > 0. There exists K = K(R, T ) < ∞ such that if where the term w I[w] − ∂ 2 xx w has been estimated as in the proof of Lemma 3.10. We infer where we have used Sobolev embedding and Lemma 3.11. Gronwall lemma yields the result.

LOCAL ERROR ESTIMATE
As pointed out in the introduction, the formalism of Lie derivatives does not seem to be easy to use in the framework of this paper, since the action of the nonlocal operator I on nonlinearities is rather involved. As a consequence, we prove an L 2 error estimate by a rather pedestrian way.
Proposition 4.1 (L 2 local error estimate). Let u 0 ∈ H 3 (R). There exists C u 0 L 2 (R) such that for all t ∈ [0, 1], . Proof. From the definition of Z t L and Remark 2.3, we have Thus, from Duhamel formula for the Fowler equation (1.4) and the Lie formula (4.1), we have: where the remainder R(t) is written as Then, from Proposition 2.1, Corollary 3.9 and Lemma 3.10, we have, for t ∈ [0, 1]: where C is a positive constant. To estimate the remainder, we decompose it as follows Let us first study the term T 1 . From Corollaries 3.3 and 3.9, we have (recall that t ∈ [0, 1]) In the same way, from Lemma 3.4 and Corollary 3.8, we control the term T 2 as From Lemma 3.2 and Corollary 3.8, For the term T 4 , write . By linearity of the evolution operator X t , we have Now from Sobolev embedding, Lemma 3.2 and Corollary 3.8, we get: Finally, since R 1 (s) = T 1 + T 2 + T 3 + T 4 then for 0 s t 1, , and by integration for s ∈ [0, t], We conclude by applying the modified fractional Gronwall Lemma 2.5.

PROOF OF THEOREM 1.2
The proof follows the same lines as in [15,Section 5]. Denote by u k = Z ∆t L k u 0 the numerical solution, and u k n = S (n−k)∆t u k , which corresponds to the exact evolution of the numerical value u k at time t k = k∆t up to time t n = n∆t. From Lemma 3.10, there exists ρ such that We prove by induction that there exists γ, ∆t 0 , c > 0 such that if 0 < ∆t ∆t 0 , for all n ∈ N with n∆t T , where C 0 = e cT u 0 H 3 . The above properties are satisfied for n = 0. Let n 1, and suppose that the induction assumption is true for 0 k n − 1. Since u n = u n n and u(t n ) = u 0 n , we estimate For k n − 2, Z ∆t L u k = u k+1 and Proposition 3.12 yields, along with the induction assumption, which is bounded by 2ρ if 0 < ∆t ∆t 0 1. Up to replacing K with max(K, 1), we obtain, for k n − 1 and n∆t T , Using Proposition 4.1, we infer for some uniform constant C. Therefore, which yields the first two estimates of the induction, provided one takes γ = CT Ke 2cT , which is uniform in n and ∆t. Finally, the last estimate of the induction stems from Corollary 3.9.

NUMERICAL EXPERIMENTS
The aim of this section is to numerically verify the Lie method convergence rate in O (∆t) for the Fowler equation (1.1). To solve the linear sub-equation (1.6), discrete Fourier transform is used and for the nonlinear sub-equation (1.7), different numerical approximations can be used. Here, we use the finite difference method. Since the discrete Fourier transform plays a key role in these schemes, we briefly review its definition, which can be found in most books. In some situation, when the mesh nodes number N is chosen to be N = 2 p for some integer p, a fast Fourier transform (FFT) algorithm is used to further decrease the computation time. In this work we will use a subroutine implemented in Matlab. In this program, the interval [0, 1] is discretized by N equidistant points, with spacing ∆x = 1/N . The spatial grid points are then given by x j = j/N , j = 0, ..., N . If u j (t) denotes the approximate solution to u(t, x j ), the discrete Fourier transform of the sequence {u j } N −1 j=0 is defined bŷ u j e −2iπjk/N , for k = 0, · · · , N − 1, and the inverse discrete Fourier transform is given by for j = 0, · · · , N − 1. Here F d denotes the discrete Fourier transform and F −d its inverse.
In what follows, the linear equation (1.6) is solved using the discrete Fourier transform and time marching is performed exactly according to where v is an average value of u in the neighbourhood of (t n , x j ).
Remark 6.1. In the case where the linear sub-equation (1.6) is solved using a finite difference scheme instead of a FFT computation, an additional stability condition is required, see [3]. Moreover, the computation time becomes very long because of the discretization of the nonlocal term which is approximated using a quadrature rule. Indeed, in [3], the Fowler equation has been discretized using finite difference method and the numerical analysis showed that this operation is computationally expensive. This observation has also motivated the use of splitting methods, in particular the implementation of split-step Fourier methods.
In order to avoid numerical reflections due to boundaries conditions and to justify the use of the FFT method, we consider initial data with compact support displayed in Figure 1 to perform numerical simulations. Since we do not know the exact solution of the Fowler equation, a classical numerical way to determine the convergence numerical order of schemes is to plot the logarithm of the error ||u 1 (T ) − u 2 (T )|| L 2 in function of the logarithm of the step time ∆t, where u 1 and u 2 are computed for time steps ∆t/2 and ∆t/4, respectively, up to the final time T . Hence, the numerical order corresponds to the slope of the curve, see Figures 2, 3. For reference, a small line of slope one (resp. two) is added in Figure 2 (resp. 3). We see that the slopes for the three initial data match well and so we can conclude that numerical simulations are consistent with the theoretical results established above.

FIGURE 2. Lie method
We also study numerical convergence of Strang splittings using initial data displayed in Figure 1. Results are plotted in Figure 3. We can see that the Strang formulation is of order two in time for smooth initial data.

FIGURE 3. Strang method
From numerical simulations, we emphasize the fact that both formulas defining a Lie operator, as well as both formulas defining a Strang operator, lead to the same results.