On the spectral vanishing viscosity method for periodic fractional conservation laws

We introduce and analyze a spectral vanishing viscosity approximation of periodic fractional conservation laws. The fractional part of these equations can be a fractional Laplacian or other non-local operators that are generators of pure jump L\'{e}vy processes. To accommodate for shock solutions, we first extend to the periodic setting the Kru\v{z}kov-Alibaud entropy formulation and prove well-posedness. Then we introduce the numerical method, which is a non-linear Fourier Galerkin method with an additional spectral viscosity term. This type of approximation was first introduced by Tadmor for pure conservation laws. We prove that this {\em non-monotone} method converges to the entropy solution of the problem, that it retains the spectral accuracy of the Fourier method, and that it diagonalizes the fractional term reducing dramatically the computational cost induced by this term. We also derive a robust $L^1$-error estimate, and provide numerical experiments for the fractional Burgers' equation.

Integro-PDEs like (1.1) typically model anomalous convection-diffusion phenomena. When µ = π λ is defined by dπ λ (z) = c λ |z| d+λ dz, c λ > 0 and λ ∈ (0, 2), (1.3) then L = −(−∆) λ/2 and the equation finds applications in e.g. over-driven detonation in gases [14] and anomalous diffusion in semiconductor growth [38]. Applications in dislocation dynamics, hydrodynamics and molecular biology can also be found, see e.g. the references in [1,18]. Many more applications can be found if asymmetric measures µ are allowed. For example, we cover the (linear) option pricing equations for all Lévy models used in mathematical finance [15,34], if dµ(z) = g(z) dπ λ (z), (1.4) for some possibly asymmetric locally Lipschitz continuous function g(·). An example is the one-dimensional (d = 1) CGMY model where g(z) = Ce −G|z| for z > 0, for positive constants C, G, M (and Y = λ). In general the non-local operator L is the generator of a pure jump Lévy process, and conversely, any Lévy process will have generator like L when (A.2) is satisfied. We refer to the books [3,15,32] for more information about Lévy processes and their many applications. The most general Lévy measures for which the results of this paper applies, are Lévy measures µ that can be decomposed as µ = µ s + µ n , (1. 5) where µ s , µ n ≥ 0, µ s is symmetric andˆ| z|>0 |z| ∧ 1 dµ n (z) < ∞. (1.6) See Section 8 for statements of results and remarks. This class possibly includes all Lévy measures, but we have so far not found a proof of this. At least it includes all the Lévy measures found in finance, see Remark 8.3, and also many singular measures like e.g. delta-measures.
It is important to note that non-linear equations like (1.1) do not admit classical solutions in general, and that shock discontinuities can develop even from regular initial conditions. This is well known for pure conservations laws (where L = 0), see e.g. [23]. For fractional conservation laws where L = −(−∆) λ/2 , it is shown in recent works that solutions are smooth for λ ∈ [1, 2) [8,18,26]. However, when λ ∈ (0, 1), the fractional diffusion is too weak to prevent shock discontinuities from forming, see [1,10,26]. In some cases however, these shocks are smoothed out over time [9]. When shocks form, weak solutions become non-unique and entropy conditions are needed to select the physically correct solution -the entropy solution. The well known Kružkov entropy solution theory for conservation laws was extended to fractional conservation laws in [1]. This extension relies on new ideas for the fractional term and is strongly influenced by the viscosity solution theory for fractional Hamilton-Jacobi-Bellman equations. Extensions of the Kružkov-Alibaud theory to general Lévy operators and even non-linear fractional terms can be found in [12,25].
In this paper we deal with a SVV (spectral vanishing viscosity) approximation of Λ-periodic entropy solutions of (1.1). The method is a Fourier Galerkin method with an additional spectral viscosity term. Because of the formation of shocks in the solutions of (1.1), it is very difficult to devise a convergent and spectrally accurate numerical approximation of this equation. This has to do with the fact that Fourier spectral methods support spurious Gibbs oscillations, and thus fails to converge strongly toward discontinuous solutions. It is well known that such methods need to be augmented by some kind of vanishing viscosity in order to achieve convergence. But the standard vanishing viscosity method is not spectrally accurate. To overcome these problems, we use the SVV approximation developed by Tadmor in [35], cf. also [11,30,33,36] and the books [4,6]. To suppress spurious oscillations without sacrificing the overall spectral accuracy of the method, Tadmor adds a modified viscosity term, which in Fourier space only affects high frequencies. There are two parameters involved in this approach, the coefficient of the viscosity term ε and the size m of the viscosity free spectrum. Spectral accuracy and convergence toward the unique, possibly discontinuous, entropy solution, then follows by imposing appropriate conditions on ǫ and m. We also like to mention another important feature of the method. In all cases, it diagonalizes the fractional term and hence reduces dramatically the computational cost induced by this term. In our rather naive implementation for the fractional Burgers' equation, the SVV method turned out to be orders of magnitude faster than a Discontinuous Galerkin approximation of the same equation where the fractional term gives full matrices.
When equation (1.1) is linear, f (u) = u, or when it is local L = 0, there is a vast literature on numerical methods and analysis, some methods and many references can be found e.g. in [4,15,23,34]. In the general case however, there is not much work on numerical methods, we only know of the papers [13,17,21]. Difference methods are introduced in [17] for equation (1.1), and in [21] for an equation similar to (1.1) from radiation hydrodynamics. In [17], the first general convergence result for monotone schemes is obtained. Finally, in [13], a Discontinuous Galerkin approximation of (1.1) is analyzed and a Kuznetsov type of theory is established and used to derive error estimates. A periodic extension of this theory will be used to find error estimates in this paper.
Throughout the paper we will use the following additional notation. A subscript p indicates Λ-periodicity in the space variables (i.e. in L ∞ p or C ∞ p ). Here Λ-periodic means 2π-periodic in each coordinate direction. As a generic constant we use C. Note that the value of C may change from line to line and expression to expression. We also need notation for high order derivatives and their norms. Let α = (α 1 , . . . , α d ) be a multi index, then Remember that α j ≥ 0, |α| = α 1 + · · · + α d , and that The rest of this paper is organized as follows. In Section 2 we introduce an entropy formulation for periodic solutions (1.1), and give a L 1 -contraction and uniqueness result. In the same section we introduce the classical vanishing viscosity approximation of (1.1), and show convergence towards (1.1) with optimal L 1 error estimate. As a corollary we get existence for (1.1). The proofs rely on the Kružkov's doubling of variables device [1,27] and Kuznetsov type of arguments [13,28], and are given in the Appendix. The SVV approximation of (1.1) is introduced in Section 3, and we show that it is spectrally accurate and that it diagonalizes the non-local operator. In sections 4-6, we assume that the measure µ is symmetric. In Section 4 we prove an energy estimate for the SVV method. Along with results from [11], this allows us to control the so-called "truncation error", the spectral projection error coming from the non-linear term. In Section 5 we prove a priori L ∞ , BV , and time regularity estimates for the SVV method, and obtain compactness. In Section 6 we prove that the SVV method converge to the classical vanishing viscosity method from Section 2. Combined with the results of that section, it follows that the SVV method converges to the entropy solution of (1.1). In the process, we also prove the optimal L 1 -rate of convergence for our SVV approximation. We solve numerically, using our SVV method, the the fractional Burgers' equation in Section 7. Finally, in Section 8 we extend the results in the previous sections to allow for asymmetric measures µ.

Entropy formulation for periodic solutions
In this section we introduce an entropy formulation for Λ-periodic solutions of the initial value problem (1.1). To this end, we write the operator L µ [·] as If r > 1, we take γ r µ = 0. The adjoint of L µ [·] takes the form We also let η, η ′ , and q denote the functions We now define the solution concept we will use in this paper.
Remark 2.1. In the entropy inequality (2.1) it is easy to see that all the terms, except possibly the L µ,r -term, are well defined and Λ-periodic in view of i). The problem with the L µ,r -term is that we integrate a Lebesgue measurable function w.r.t. a Radon measure µ. But the term is still well defined because the integrand of L µ,r [u] is measurable w.r.t. the product measure dµ(z)dxdt. This is true because the integrand is the dµ(z)dxdt-a.e. limit of continuous functions, a fact which readily follows from the fact that u is the dxdt-a.e. limit of smooth functions. By i), (A.3), and Fubini, we then find that L µ,r [u] ∈ C([0, T ]; L ∞ p (R d )). We now state the following central result: Theorem 2.2. (L 1 -contraction) Let u and v be two entropy solutions of the initial value problem (1.1) with initial data u 0 and v 0 . Then, for a.e. t ∈ (0, T ), The proof will be given in Appendix A. Uniqueness for periodic entropy solutions of (1.1) immediately follows by setting u 0 = v 0 .

Corollary 2.3. (Uniqueness)
There is at most one entropy solution of (1.1).
We now consider the vanishing viscosity approximation of (1.1), In this paper we always assume that this problem admits a unique classical solution u ǫ . This is of course true, but a proof lays outside the scope of this paper. Remark 2.6 below provides some ideas on how to prove this result. We now give an estimate on the rate of convergence of u ε toward the entropy solution u of (1.1).
Theorem 2.4 (Convergence rate I). Let u be the periodic entropy solution of (1.1), and u ǫ be a smooth solution of (2.3). Then, The proof is given in Appendix B. This result generalizes to periodic fractional conservation laws Kuznetsov's well known result for scalar conservation laws [28]. As a by-product of the well-posedness of (2.3) and Theorem 2.4, we have the existence of entropy solutions of (1.1).
Remark 2.6. Uniqueness of solutions of (2.3) can be proved using an entropy formulation (see the start of Appendix B) and a standard adaptation of the proof of Theorem 2.2 incorporating ideas of Carrillo [7] to handle the Laplace term. Existence of an entropy solution can be proven e.g. by appropriately modifying our spectral approximation, compactness, and convergence analysis, see the following sections. The solution of (2.3) will also be smooth. To see this, note that the principal term in (2.3) is the ǫ∆-term while the L µ -term is of lower order, and hence regularity proofs for viscous conservation laws ((2.3) with µ ≡ 0 and ε > 0) should still work after some modifications. We refer to e.g. [31] for regularity of viscous conservation laws, and note that the modifications typically consist of using interpolation inequalities for the L µ -term, see e.g. Lemma 2.2.1 in [22].

The spectral vanishing viscosity method
We introduce a Fourier spectral method for the d-periodic initial value problem (1.1). The approximate solutions will be N -trigonometric polynomials, which solve the semi-discrete spectral vanishing viscosity (SVV) approximation where the Fourier projection P N is defined as The (spectral) vanishing viscosity term has the following three ingredients: Combined with one's favorite ODE solver (e.g. Euler, Runge-Kutta, etc.), (3.1) and (3.2) give a fully discrete numerical approximation method for (1.1). With left-hand sides set to zero (µ ≡ 0 and ε N = 0), (3.1) becomes the standard Fourier approximation of (1.1). It is well known that this approximation is spectrally accurate but, as opposed to the equation, it lacks entropy dissipation. The approximation supports spurious Gibbs oscillations which prevent strong convergence toward solutions containing shock discontinuities. If the L µ -term is present in the equations, shock solutions are still possible in some situations [2], and the problem of the Gibbs oscillations remains. In order to suppress such oscillations without sacrificing the overall spectral accuracy of the method, we have followed Tadmor [35] and added a vanishing spectral viscosity term to the scheme, An important feature of Fourier method (3.1) is that it diagonalizes, and hence localizes, the non-local operator L µ [·]! This leads to dramatically reduced computational cost for this term. Indeed, Furthermore, when the measure µ is symmetric, the weights (3.4) are all real and non-positive. This follows since the imaginary part of the integrand is odd and the real part is even and non-positive (e iξ·z = cos(ξ · z) + i sin(ξ · z)). Finally, we stress that the approximation of the non-local operator (1.2) is spectrally accurate since, by Taylor's formula, and note that To conclude this section, we recall that by Lemma 3.1 and Corollary 3.2 of [11], the spectral vanishing viscosity term is an L p -bounded perturbation of the standard vanishing viscosity ǫ N ∆u N :

Spectrally small truncation error for symmetric µ
In this section we assume that the measure µ is symmetric, cf. (3.5). In the SVV approximation (3.1), the convection term We will now show that this error is spectrally small due to the presence of the spectral vanishing viscosity term.
Let us start by noting that a straightforward estimate leads to for all multi-indices α, β. Note that there is no divergence in this estimate, so ∂ α x f is a vector. By Theorem 7.1 in [11], there is a constant K s such that ΩN ) and Ω N = {u : |u| ≤ u N L ∞ (Λ) }. This inequality is a type of Gagliardo-Nirenberg-Moser estimate, and similar results can be found in page 22 in Taylor [37]. By these two inequalities we can conclude that, for all 0 ≤ r ≤ s, Inequality (4.2) states that the r-derivative of the truncation error decays as rapidly as the s-smoothness of u N permits. Of course the s-derivatives of an arbitrary N -trigonometric polynomial u N may grow as fast as N s , in which case nothing is gained from (4.2). However, if u N is solves our VVS approximation (3.1), we can have the better bound ǫ −s N in L 2 . This will be a consequence of the following energy estimate: Remember that in this section µ is symmetric and hence G µ is real and nonpositive. Now if (A. 8) |f | C s < ∞ for sufficiently large s, cf. (4.7) below, and Taking into account (4.2), we then find that We can now turn these inequalities into spectral decay estimates in the uniform norm using the Sobolev inequality (cf. Theorem 6, Chapter 5, in [20]) For example, inequality (4.5) becomes Note that the polynomial decay rate in (4.6) can be made as large as the C ssmoothness of f (·) permits. Taking r = 2, we can find the following result.
The smoothness requirement (4.7) will be sufficient for all the estimates derived throughout the paper.
Proof of Theorem 4.1. For sake of brevity, we will write · instead of · L 2 (Λ) . With (3.6) in mind, we rewrite the SVV approximation (3.1) in the two equivalent forms and hence spatial integration of (4.10) against u N yields Using (3.8) with p = 2 for the first term on the right and (4.2) with (r, s) = (0, 1) for the second term, we find that Hence (4.3) follows for s = 0 since by (A.7), . The general case follows by induction on s. Spatial integration of (4.9) against ∂ 2α x u N for some multi-index α yields 1 2 (4.11) After having used (3.8) and Young's inequality to bound the first and second term on the right hand side, we find that Now we sum over all |α| = s to find that (4.12) By (4.1) and (4.2), and hence by inequality (4.12) we see that where the last inequality holds for N big enough. By (A.7) and integration in time, we then find that (4.14) At this point (4.3) follows by the induction assumption on s since The proof is now complete.

A priori estimates and compactness
In this section we prove uniform Proof. For sake of brevity, we write just · ∞ instead of · L ∞ (Λ) . First we note that, for any smooth convex function η(·) with derivative η ′ (·), we have that This is a consequence of the inequality η ′ (b)(a − b) ≤ η(a) − η(b) which holds for all smooth convex functions η(·). Moreover, To see this note that since u N is smooth and periodic. By Fubini we then find that By Λ-periodicity of u N , (5.2) now follows sincê Λ η(u N (x + z)) dx =ˆΛ η(u N (x)) dx for every z, and Let us now integrate (4.10) against the function p u p−1 N (with p even), and use (5.1) and (5.2) to get rid of the non-local operator L µ [·]. We then find that which by the Hölder inequality (with p and q = p p−1 ) is less than or equal to Since φ p−1 L p p−1 = φ p−1 L p , we may divide both sides by p u N (·, t) p−1 L p (Λ) and send p → ∞ to discover that By (4.8), (3.8), the definitions of B s , K s and c N , and (A.7), it follows that Letting y(t) = e −cN t u N (·, t) ∞ , and multiplying by the integrating factor e −cN t , we find that dy dt Going back to u N (·, t) ∞ , we can conclude that where the last factor is bounded for t ≤ C ln N for some C.
We also have the following result: If we integrate this expression against sgn ̺ (∂ i u N ), where sgn ̺ (·) is a smooth approximation of the sign function, we can get rid of the non-local operator L µ [·] as in the proof of Lemma 5.1. If we also use (3.8) with p = 1 and take the limit as ̺ → 0, a standard computations reveal that d dt , we integrate this inequality in time to see that But by (4.4), and the proof is complete.
Proof. Let u ǫ N (·, t) = u N (·, t) * ω ǫ (·) for an approximate unit ω ǫ (cf. the proof of Theorem 2.2). By the triangle inequality we see that 3) The first and the third term on the right-hand side of (5.3) are bounded by ǫ|u| BV : Let us estimate the second term. By Taylor's formula with integral remainder, We now derive a bound for ∂ t u N L 1 (and hence also for ∂ t u ε N L 1 ) by using the SVV approximation (3.1) itself. To this end, we take the convolution product of both sides of (4.9) with ω ǫ to obtain By the triangle inequality and Young's inequality for convolutions, Therefore, by the regularity of f and u N ((A.8), Lemmas 5.1-5.2) and (4.8), we find that For the term containing the non-local operator we write The second term on the right hand side of the inequality above is easily seen to be bounded by C u N (·, t) L 1 , while Taylor's formula with integral reminder and integration by parts reveals that the first term is bounded bŷ For the Laplace term we have and finally, using Young's inequality for convolutions and (3.8), To sum up we have and inequality (5.3) and the above estimates then implies that Take ε = |t 1 − t 2 | and the proof is complete.

Convergence and error estimate
The solution v ǫN of the vanishing viscosity method (2.3) converges to the unique entropy solution u of (1.1), and by Theorem 2.4, In this section we prove a similar error estimate between v ǫN and the SVV approximation u N .
A direct consequence of Theorems 2.4 and 6.1, is the following convergence and error estimate for the SVV method.
Proof of Theorem 6.1.
Since v ǫN is smooth, we can subtract equation (2.3) from equation (3.1) to obtain As explained in the proof of Lemma 5.2, we can integrate such an inequality against (a smooth approximation of) sgn(u N − v ǫN ), to find that (after going to the limit) By By (4.5), . The proof is now complete.

An application: the fractional Burgers' equation
In this section we apply the results of the previous sections to numerically solve the fractional (or fractal) Burgers' equation in R d , where the fractional Laplacian term −(−∆) λ/2 u N = L π λ [u N ] and π λ has been defined in (1.3). In this setting expression (3.4) becomes , cf. [19]. We have the following result: Proposition 7.1.
. The proof is given at the end of this section. In the above result and in the following, dS y will denote the surface area measure of the unit sphere |y| = 1. Expression (7.2) is the "Fourier symbol" of the fractional Laplace operator in our periodic setting. When λ ∈ (0, 1), the integral Θ λ =´∞ 0 x −λ sin x dx is a generalized Fresnel integral [29] with value When λ = 1, Θ λ is a Dirichlet integral [24] and has value π 2 . For λ ∈ (1, 2), the integral Θ λ has to be evaluated numerically since explicit formulas are not available.
Simple energy estimates can then be used to show that the solutions of (7.1) belong to H λ/2 (Λ), which is more regularity than what can be expected for general solutions of the pure Burgers' equation (µ = 0). We now use the SVV method (3.1) to work out some approximate solutions of the fractional Burgers' equation (7.1) with d = 1. Hence f (u) = u 2 /2 and µ = π λ in (3.1). We multiply both sides of (3.1) by e −iξ x , and integrate over (0, 2π) to obtain the following system of ODEs where the Fourier coefficientsQ ξ satisfy the assumptions listed in Section 3 and are chosen as in [30] (they vary continuously between zero and one). In our simulations we have used a fourth order Runge-Kutta solver for (7.3).
The results of our numerical simulations can be found in Figure 1 and 2. The results in Figure 1 confirm the convergence of our SVV approximation (7.3) for all all values of λ ∈ (0, 2). In Figure 2 we have solved the (7.3) with ǫ N = 0 (no spectral vanishing viscosity). For λ > 1, convergence continues to hold, while for λ < 1, convergence fails and spurious Gibbs oscillations appear. This is consistent with the theoretical results for fractional conservation laws [2,18]: These equations admit smooth solutions for λ > 1 (the strong diffusion case), while shock discontinuities may appear for λ < 1 (the weak diffusion case).
Proof of Proposition 7.1. Let us prove the case d = 1 first. By Euler's formula, e iξz = cos(ξz) + i sin(ξz), we find that Taylor expansions show that these integrals are finite. In fact, the sin-integral is zero since the its integrand is odd. Integration by parts then leads tô Now we consider the integral over |z| > r. Again the imaginary part (the sine part) is zero, and a computation like the one we performed above reveals that Note the + sign of the cosine-term! We add these two equations and find that The integral´∞ 0 z −λ sin(ξz) dz is finite and positive for all λ ∈ (0, 2) (cf. [16] for details). Whenever ξ > 0, we can use the change of variable ξz → x to deduce that When ξ < 0, we use the relation sin(−ξx) = − sin(ξx) to obtain G π λ (ξ) = − 2c λ λ |ξ| λˆ∞ 0 sin x x λ dx, and the conclusion for d = 1 follows. When d > 1 we use polar coordinates x = ry for r > 0 and |y| = 1, and we find thatˆ| Proceeding as in the d = 1 case for the r-integral with y fixed, we find that By symmetry, the value of the y-integral is the same for any ξ. Therefore, The proof for the case d > 1 is now complete.

Extension to asymmetric measures µ
In this section we show how to modify the arguments of the previous sections to obtain results for a large class of non-symmetric measures µ including all the Lévy measures used in finance. A careful look at the previous arguments shows that symmetry of µ is used for the sole purpose of having a sign of the fractional term in the energy inequality (see (4.14)) in order to prove Theorems 4.1 and 4.2. This fractional term is and it is non-positive when µ is symmetric. In the general case the sign of the fractional term (8.1) is unknown, but everything still works if we assume that µ = µ s + µ n , for µ s , µ n satisfying (1.6) (i.e. we assume (1.5) and (1.6)). Note that in this case, we may split the weights in (3.4) into their symmetric and non-symmetric parts, where G µs (ξ) is again real and non-positive, and by (1.6), The main result of this section is the following: Note that µ s , µ n ≥ 0, µ s is symmetric, and that µ n satisfies the integrability condition in (1.6) if g is locally Lipschitz: Let g n (z) = g(z) − g(z) ∧ g(−z) and note that g n (0) = 0, hence g n (z) = |g n (z) − g n (0)| ≤ C |z| for |z| < 1.
We now show how to modify the proof of Theorem 4.1 to prove Theorem 8.2.
Proof of Theorem 8.2. Once again we use the shorthand · instead of · L 2 (Λ) , and rewrite the SVV approximation (3.1) as in (4.9) and (4.10). Note that (5.1) and (5.2) holds for general measures µ, so we find that Hence, spatial integration of (4.10) against u N yields and the conclusion in the case s = 0 follows exactly as in the first part of the proof of Theorem 4.1. Now let s > 0, and note that by (8.2) and Young's inequality, If we take this into account and perform spatial integration of (4.9) against ∂ 2α x u N for some multi-index α, we find the following modified version of (4.11), As in the proof of Theorem 4.1, we now use (3.8) and Young's inequality to bound the first and second term on the right hand side. The result is that Now we sum over all |α| = s to find that Thanks to (4.1) and (4.2), and hence where the last inequality holds for N big enough.
Finally, the theorem follows from renaming t 2 and using part iii) in Definition 2.1 to send t 1 → 0. The results then follows by setting ρ = √ ǫ and φ = χ µ as in (A.5), and conclude as in the proof of Theorem 2.2: Sending µ → 0, setting t 2 = t, and using part iii) in Definition 2.1 to send t 1 → 0.