Fourier Spectral Methods for Some Linear Stochastic Space-Fractional Partial Differential Equations

Abstract: Fourier spectral methods for solving some linear stochastic space-fractional partial differential equations perturbed by space-time white noises in the one-dimensional case are introduced and analysed. The space-fractional derivative is defined by using the eigenvalues and eigenfunctions of the Laplacian subject to some boundary conditions. We approximate the space-time white noise by using piecewise constant functions and obtain the approximated stochastic space-fractional partial differential equations. The approximated stochastic space-fractional partial differential equations are then solved by using Fourier spectral methods. Error estimates in the L2-norm are obtained, and numerical examples are given.

Let us here consider two examples, which apply the fractional Laplacian in the physical models.The first example is the surface quasi-geostrophic (SQG) equation, where κ ≥ 0 and α > 0, θ = θ(x 1 , x 2 , t) denotes the potential temperature, u = (u 1 , u 2 ) is the velocity field determined by θ.When κ > 0, the SQG equation takes into account the dissipation generated by a fractional Laplacian.The SQG equation with κ > 0 and α = 1/2 arises in geophysical studies of strongly-rotating fluids.For the dissipative SQG equation, α = 1/2 appears to be a critical index.In the subcritical case when α > 1/2, the dissipation is sufficient to control the nonlinearity, and the global regularity is a consequence of a global a priori bound.In the critical case when α = 1/2, the global regularity issue is more delicate.There are few theoretical results for the supercritical case α < 1/2 in the literature [7].
The second example is about the wave propagation in complex solids, especially viscoelastic materials (for example, polymers) [8].In this case, the relaxation function has the form k(t) = ct −ν , 0 < ν < 1, c ∈ R, instead of the exponential form known in the standard models.This polynomial relaxation is due to the non-uniformity of the material.The far field is then described by a Burgers equation with the leading operator (−∆) 1+ν 2 instead of the Laplacian: This equation also describes the far-field evolution of acoustic waves propagating in a gas-filled tube with a boundary layer.
Frequently, the initial value or the coefficients of the equation are random; therefore, it is natural to consider the stochastic space-fractional partial differential equations.The existence, uniqueness and regularities of the solutions of stochastic space-fractional partial differential equations have been extensively studied; see, for example, [3,4,9,10].In this work, we will focus on the case 1/2 < α ≤ 1, since the existence, uniqueness and regularity of the solution in this case is well understood in the literature; see [11] (Theorem 1.3).However, the numerical methods for solving space-fractional stochastic partial differential equations are quite restricted even for the case 1/2 < α ≤ 1. Debbi and Dozzi [11] introduced a discretization of the fractional Laplacian and used it to obtain an approximation scheme for the fractional heat equation perturbed by a multiplicative cylindrical white noise.As far as we know, [11] is the only existing paper in the literature that deals with this kind of numerical approach for such a problem.In this work, we will use the ideas developed in [12] to consider the numerical methods for solving stochastic space-fractional partial differential equations; see also [13][14][15][16].We first approximate the space-time white noise by using piecewise constant functions and then obtain the approximate solution û(t) of the exact solution u(t).Finally, we provide error estimates in the L 2 -norm for u(t) − û(t).
For the deterministic space-fractional partial differential equations, many numerical methods are available in the literature.There are two ways to define the fractional Laplacian.One way of defining (−∆) α v, 1/2 < α ≤ 1 is by using the eigenvalues and eigenfunctions of the Laplacian −∆ subject to the boundary conditions as in (4).Another way of defining (−∆) α v, 1/2 < α ≤ 1 is by using the integral for the function ṽ, where ṽ is defined on the whole real line R and is the extension function of v: More precisely, for ṽ(x), x ∈ R, we define: where C α is a positive constant depending on α.We then define [17], where F and F −1 denote the Fourier and inverse Fourier transforms, respectively.For v(x), x ∈ (0, 1), we define the fractional Laplacian by: It is easy to show that for some suitable functions w(x), x ∈ R [18], Hence, for the function v(x) defined on the bounded interval (0, 1), we have: which is also called the Riesz fractional derivative.We note that Definitions (4) and ( 5) are not equivalent [17].For the deterministic space-fractional partial differential equations where the space-fractional derivative is defined by (5), or the Riemann-Liouville space-fractional derivative, or the Caputo space-fractional derivative, many numerical methods are available, for example, finite difference methods [18][19][20][21][22][23][24][25][26][27][28][29][30], finite element methods [14,[31][32][33][34][35][36][37][38][39][40] and spectral methods [41,42].For the deterministic space-fractional partial differential equations where the space-fractional derivative is defined by (4), some numerical methods are also available, for example the matrix transfer method (MTT) [21,22,43] and the Fourier spectral method [44].In this work, we will use Fourier spectral methods to solve the approximated stochastic space-fractional partial differential equations.The main advantage of this approach is that it gives a full diagonal representation of the fractional operator, being able to achieve spectral convergence regardless of the fractional power in the problem.Let 0 = x 0 < x 1 < x 2 < • • • < x J = 1 be the space partition of (0, 1) and h the space step size.Let 0 = t 0 < t 1 < t 2 < • • • < t N = T be the time partition of (0, T) and k the time step size.To find the approximate solution of (1)-(3), we first approximate the space-time white noise ∂t∂x by using a piecewise constant function defined by, with n = 1, 2, 3, ..., N, j = 1, 2, ..., J [12], where η nj ∈ N (0, 1) is an independently and identically distributed random variable and: Hence: We also note that [12]: The solution u(t, x) of ( 1)-( 3) can then be approximated by û(t, x), which solves the following: Note that is a function in L 2 ((0, T) × (0, 1)), and therefore, we can solve ( 8)-(10) by using any appropriate numerical method for deterministic space-fractional partial differential equations.In Theorem 2, we prove that, if 1/2 < α ≤ 1, then: Let us now introduce the Fourier spectral method for solving (8)- (10).Let J be a positive integer, and denote: S J = span{e 1 , e 2 , . . ., e J } Define by P J : H → S J the projection from H to S J , The Fourier spectral method for solving (8)-( 10) is to find ûJ (t) ∈ S J , such that: In Theorem 4, we prove that: Combining Theorem 2 with Theorem 4, we have: The paper is organized as follows.In Section 2, we consider the approximation of space-time white noise.In Section 3, we consider a Fourier spectral method for deterministic space-fractional partial differential equations, and the error estimates are proven.In Section 4, we provider a numerical example.
Note that t ≥ t j , we then have, by using (63), For I 2 , we have: where I 21 can be estimated as follows: which implies, by (64): For I 22 , we have: Moreover, applying (66) and taking into account | cos(nπx k+1 ) − cos nπx k | ≤ nπh, we derive: For I I we have, with η kj = N (0, 1), Similarly, we can estimate I I I. Together, these estimates complete the proof of Theorem 3.
Similarly, the solution of ( 34)-( 36) has the form of: Theorem 4. Assume that û and ûJ are the solutions of (33) and (37), respectively.Let 0 ≤ r < 1/2, and assume that u 0 ∈ H r 0 (0, 1).Then, there exists a positive constant C, such that: In particular, we have, with r = 0, (39) To prove Theorem 4, we need the following smoothing property for the solution operator E α (t).
1. Let s > 0. We have: for some constants C and δ which depend on s and α. 2. Let P J : H → S J be defined by (12), then: Proof.Note that A is a positive definite operator with eigenvalues 0 < λ 1 < λ 2 < λ 3 < . . . .For any function g(•), we have: which is (40).To show (41), we note that: The proof of Lemma 5 is complete.

Numerical Simulations
In this section, we will present the computational issues for solving the following stochastic space-fractional parabolic partial differential equations by using the spectral method developed in the previous section, with 1/2 < α ≤ 1, where (−∆) α is the fractional Laplacian defined by using the eigenvalues and eigenfunctions of the Laplacian operator −∆ subject to some boundary conditions.Here, f : R → R is a smooth function, and > 0 denotes the diffusion coefficient.In our numerical example, we will use the discrete sine transform MATLAB functions dst and idst.We also include the nonlinear term f , although the error estimates in the previous sections are only proven for f = 0.In our future work, we will consider the error estimates for solving the nonlinear stochastic space-fractional partial differential equations with multiplicative noise by using the spectral method.Let x 0 < x 1 < • • • < x J = 1 be a space partition of [0, 1] and ∆x = h be the space step size.Let 0 = t 0 < t 1 < • • • < t N = T be the time partition of [0, T] and ∆t = k the time step size.The space-time noise ∂ 2 W(t,x) ∂t∂x is approximated by using piecewise constant function ∂t∂x , where: For convenience, we will denote Ĝ(t, x) = ∂ 2 Ŵ(t,x) ∂t∂x below.Equations ( 43)-(45) can then be approximated by the following, with 1/2 < α ≤ 1, . Then, A has eigenvalues λ j and eigenfunctions e j where: That is: Ae j = λ j e j , j ∈ Z + Equations ( 47)-( 49) can further be written as the following abstract form: find û(t) ∈ H 1 0 (0, 1) ∩ H 2 (0, 1), such that: Let S J−1 := span{e 1 , e 2 , . . ., e J−1 }.The spectral method for solving (47)-( 49) is to find u J−1 (t) ∈ S J−1 , such that, with 0 < t ≤ T, where P J−1 : H → S J−1 is the orthogonal projection operator defined by: where A J−1 = P J−1 A : S J−1 → S J−1 and (•, •) denotes the inner product in H = L 2 (0, 1).We remark that we use S J−1 (not S J ), since we will apply the MATLAB functions dst and idst in our numerical algorithms below.The semi-implicit Euler method for solving (47)-( 49) is to find It is easy to see that the Fourier coefficients ũj,n satisfy, with j = 1, 2, . . ., J − 1, ũj,0 = (P J−1 u 0 , e j ) where: Here, ũj,n , Gj,n , fj ( ûJ,n ) denote the Fourier coefficients of ûJ−1,n , Ĝ(t n ) and f ( ûJ−1,n ), respectively.We may use the following steps to describe how to solve (47)-(49) numerically by using the spectral method: Step 1: Given initial value û0 (x) and f , we get the approximation u J−1,0 (x) = P J−1 u 0 ≈ u 0 and Step 2: Find the Fourier coefficients ũj,0 and fj ( ûJ−1,0 ) by: and: Here, ), and Ŵ is generated by: Step 3: Find the Fourier coefficients ũj,1 , j = 1, 2, . . ., J by: where ./denotes the element-wise division and: Step 4: Find the Fourier coefficients ũj,2 , j = 1, 2, . . ., J − 1 by: ûj,2 = (1 + ∆tλ j ) −1 ũj,1 + ∆t fj ( ûJ−1,1 ) + ∆t Gj,1 where: ), and Ŵ is defined in (59).
Let us now introduce the MATLAB functions to solve our problem.Let u 0 denote the initial value vector, that is u 0 = [u 0 (x 1 ), u 0 (x 2 ), . . ., u 0 (x J−1 )].Let u denote the approximate solution vector at time T, that is u = [u(x 1 , T), u(x 2 , T), . . ., u(x J−1 , T)].We may use the following MATLAB function to get the approximate solution u at T for any function f .Here, we choose We can obtain the approximate solution u at time T at the different x k , k = 1, 2, . . ., J − 1 by the following MATLAB function.
In Figure 2, we plot an approximation sample path of u(t, x) with J = 2 4 and N = 2 6 on 0 ≤ t ≤ 1 and 0 ≤ x ≤ 1.In Figure 3, we consider the convergence rate against the different time steps.Choose the fixed J = 64; we then consider the different time steps.The reference solution is obtained by using the time step ∆tre f = T/Nre f with Nre f = 10 4 .Let kappa = [20,50,100,150,200,250, 300]; we will consider the approximate solutions with the different time steps ∆t i = ∆tre f * kappa(i), i = 1, 2, . . ., 7. In our experiment, for saving the computation time, we will consider the error estimates ûN (t n ) − u(t n ) L 2 (Ω,H) at time t n .We hope to observe the same convergence order as in Theorem 6.
To do this, we consider M = 100 simulations.For each simulation ω m , m = 1, 2, . . ., M, we compute ûN (t n ) ≈ û(t n ) at time t n = 1 by using the different time steps.We then compute the following L 2 norm of the error at t n = 1 for the simulations ω m , m = 1, 2, . . ., M, where the reference (or "true ") solution uref(t n , ω m ) is approximated by the time step ∆tre f = T/Nre f .We then average (∆t i , ω m ) with respect to ω m to obtain the following approximation of ûN (t n )) − uref(t n ) L 2 (Ω,H r ) with respect to the different time step ∆t i , Since the convergence rate with respect to the time step is O(∆t 1/2 ), i.e., S(∆t i ) ≈ ∆t 1/2 i this implies that: log(S(∆t i )) ≈ 1/2 log(∆t i ), i = 1, 2, . . ., 7 In Figure 3, we plot the points log(∆t i ), log(S(∆t i )) , i = 1, 2, . . ., 7, and we see that the points are parallel to the reference line, which has the slope 1/2, as we expected in our theoretical results.
In Table 2, we list the error S(∆t i ) against the different time steps ∆t i .In Figure 4, we plot the L 2 error S(∆t) against the different J where the L 2 errors are approximated by using M = 100 simulations.We indeed observe the convergence with respect to the different J.

Conclusions
In this work, we present a Fourier spectral method for solving space-fractional partial differential equations.The space-time white noise is approximated by using piecewise constant functions.For the linear problem, we obtain the exact error estimates in the L 2 -norm and find the relations between the convergence order and the fractional power α, 1/2 < α ≤ 1.For the nonlinear problem, we introduce the numerical algorithm and the MATLAB code for solving it based on the discrete sine transform and inverse discrete sine transform MATLAB functions dst.m and idst.m.The MATLAB code in this paper can be easily modified to solve other nonlinear stochastic fractional partial differential equations with Dirichlet boundary conditions.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In this Appendix, we shall provide the proof of Theorem 2. To do this, we need the following lemmas.
We also need the following isometry property for space-time white noise W(s, y); see, e.g., [1].

Figure 4 .
Figure 4.A plot of the error at T = 1 against the J.