Sharp convergence rates of time discretization for stochastic time-fractional PDEs subject to additive space-time white noise

Sharp convergence rates of time discretization for stochastic time-fractional PDEs subject to additive space-time white noise

Max Gunzburger Department of Scientific Computing, Florida State University, Tallahassee, FL 32306, USA. (gunzburg@fsu.edu)    Buyang Li Department of Applied Mathematics, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong. (buyang.li@polyu.edu.hk)    Jilu Wang Department of Scientific Computing, Florida State University, Tallahassee, FL 32306, USA.(jwang13@fsu.edu)
Abstract

The stochastic time-fractional equation with space-time white noise is discretized in time by a backward-Euler convolution quadrature for which the sharp-order error estimate

is established for , where denotes the spatial dimension, the approximate solution at the time step, and the expectation operator. In particular, the result indicates optimal convergence rates of numerical solutions for both stochastic subdiffusion and diffusion-wave problems in one spatial dimension. Numerical examples are presented to illustrate the theoretical analysis.
Keywords: stochastic partial differential equation, time-fractional derivative, space-time white noise, error estimates

1 Introduction

We are interested in the convergence of numerical methods for solving the stochastic time-fractional partial differential equation (PDE) problem

(1.1)

where , , denotes a bounded region with Lipschitz boundary , a given deterministic source function, given deterministic initial condition, and a space-time white noise, i.e., the time derivative of a cylindrical Wiener process in . The underlying probability sample space for the stochastic noise is denoted by . The operator denotes the Laplacian, defined on the domain

and denotes the left-sided Caputo fractional time derivative of order , defined by (c.f. [10, pp. 91])

(1.2)

where denotes the Euler gamma function.

Problem (1.1) arises naturally by considering the heat equation in a material with thermal memory, subject to stochastic noise [3, 11, 16]. For the model (1.1), both the fractional time derivative and the stochastic process forcing result in solution having low regularity. Hence, the numerical approximation of such problems and the corresponding numerical analysis are very challenging. By defining for and using the identity

(1.3)

applying to (1.1) yields another formulation of (1.1):

(1.4)

where the case can be verified directly from (1.1). For the sake of clarity, we focus on only one of the equivalent problems, namely (1.1).

The solution of (1.1) can be decomposed into the solution of the deterministic problem

(1.5)

plus the solution of the stochastic problem

(1.6)

The stability and convergence of numerical solutions of (1.5) have been widely studied [2, 4, 14, 15, 17]. For example, if is smooth in time then numerical methods of up to order are available for approximating the solution of (1.5) and its equivalent formulations [8, 9, 14, 15]. In particular, the convolution quadrature generated by the backward Euler method yields a first-order convergence rate for solving (1.5).

In this work, we focus on numerical approximation of the stochastic time-fractional PDE (1.6) with additive space-time white noise based on the convolution quadrature generated by the backward Euler method. In the case and , rigorous error estimates for numerical solutions of (1.6) are carried out in [11] for the case of additive Gaussian noise in the general -Wiener process setting. For a space-time white noise, an almost optimal-order convergence rate for the time discretization error

is proved [11, Remark 4.7, with ]) for arbitrarily small , where and denote the mild solution and numerical solution of (1.6) at time , respectively. We are not aware of any rigorous numerical analyses in the case . In the special case , error estimates for time discretization of the stochastic PDE (1.6) are proved in [7] and [6, 18] for Rothe’s method and the backward Euler method, respectively. In particular, the sharp convergence rate in time is proved for the expectation of the norm time discretization error [1].

The aim of this paper is to prove, for general -dimensional domains, the sharp-order convergence rate estimate

(1.7)

for time discretization of the stochastic PDE (1.6). This estimate is achieved via a more delicate analysis of the resolvent operator by using its Laplace transform representation. Our result covers both subdiffusion and diffusion-wave cases in one-dimensional spatial domains and, for the subdiffusion case, multi-dimensional domains.

The rest of the paper is organized as follows. In Section 2, we present the backward-Euler convolution quadrature scheme we use to determine approximate solutions of the stochastic time-fractional PDE (1.6) and then state our main theoretical results. In Section 3, we derive an integral representation of the numerical solution for which we prove sharp convergence rate results for the approximate solution. Numerical results are given in Section 4 to illustrate the theoretical analyses.

Throughout this paper, we denote by , with/without a subscript, a generic constant independent of and which could be different at different occurrences.

2 The main results

In this section, we describe the time discretization scheme we use for determining approximate solutions of the stochastic time-fractional PDE (1.1) and state our main results about the convergence rate of the numerical solutions.

2.1 Mild solution of the stochastic PDE

Let , , denote the -norm normalized eigenfunctions of the Laplacian operator corresponding to the eigenvalues , arranged in nondecreasing order. The cylindrical Wiener process on can be represented as (cf. [5, Proposition 4.7, with and denoting some negative-order Sobolev space])

(2.1)

with independent one-dimensional Wiener processes , .

In the case , the solution of the deterministic problem (1.5) can be expressed by (via Laplace transform, cf. [14, (3.11) and line 4 of page 12])

(2.2)

where the operator is given by

(2.3)

with integration over a contour on the complex plane,

(2.4)

The angle above can be any angle such that so that, for all to the right of in the complex plane, with .

Correspondingly, the mild solution of (1.6) is define as (cf. [16] and [11, Proposition 2.7])

(2.5)
(2.6)

This mild solution is well defined in ; see the Appendix A.

2.2 Convolution quadrature

Let denote a uniform partition of the interval with a time step size , and let . Under the zero initial condition, the Caputo fractional time derivative can be discretized by the backward-Euler convolution quadrature (also known as Grünwald-Letnikov approximation, cf. [4]) defined by

(2.7)

where , , are the coefficients in the power series expansion

(2.8)

Here, is the characteristic function of the backward-Euler method and we set

(2.9)

For any sequence , we denote the generating function of the sequence by

(2.10)

that is an -valued analytic function in the unit disk and the limit

exists in . Then, we have

(2.11)

2.3 Time-stepping scheme and main theorem

With the notations introduced in the last subsection, we discretize the fractional-order derivative in (1.6) by using convolution quadrature in time to obtain

(2.12)

Equivalently, can be expressed as

(2.13)

The main result of this paper is the following theorem.

Theorem 2.1.

Let with . Then, for each the numerical solution given by (2.12) is well defined in and converges to the mild solution with sharp order of convergence, i.e., we have

(2.14)

where denotes the expectation operator and the constant is independent of .

3 Proof of Theorem 2.1

3.1 The numerical solution in

In this subsection, we show that the numerical solution is well defined in . To this end, we use the following estimate for the eigenvalues of the Laplacian operator. For the simplicity of notations, we denote by and the inner product and norm of , respectively.

Lemma 3.1 ([12, 13]).

Let denote a bounded domain in , . Suppose denotes the eigenvalue of the Dirichlet boundary problem for the Laplacian operator in . With denoting the volume of , we have that

(3.1)

for all , where and denotes the volume of the unit -ball.

Clearly, if for , then

(3.2)

In view of (2.13), we only need to prove

(3.3)

In fact, we have

Hence, for a fixed time step ,

is a Cauchy sequence in . Consequently, (3.3) is proved. In view of (2.13) and (3.2)-(3.3), the numerical solution is well defined in .

3.2 A technical lemma

To prove the error estimate in Theorem 2.1, we need the following technical lemma.

Lemma 3.2.
(3.4)
(3.5)
Proof.

Clearly, Lemma 3.1 implies .

First, if ,

(3.6)

Second, if , by setting to be the largest integer that does not exceed , we have

It is easy to see that and

where the equality follows by changing the variable . This proves (3.4) in the case .

Finally, for the point in the complex plane, we have and . By looking at the triangle with three vertices , , and with interior angles , , and at the three vertices, respectively, we have

If , then would be the length of the longest side of the triangle, i.e.,

which immediately implies

If , then the angle condition implies . Hence, we have

which immediately implies

In either case, we have (3.5). This completes the proof of Lemma 3.2. ∎

3.3 Solution representations

In this subsection, we derive a representation of the semidiscrete solution by means of the discrete analogue of the Laplace transform and generating function.

Let denote the truncated piece of the contour defined by

(3.7)

For , let denote the segment of a vertical line defined by

(3.8)

The following technical lemma will be used in this and next subsections.

Lemma 3.3.

Let and be given, and let be fixed, with defined in (2.9). Then, both and are analytic with respect to in the region enclosed by

Furthermore, we have the following estimates:

(3.9)
(3.10)
(3.11)
(3.12)

where the constants , , and are independent of and .

Proof.

Clearly, (3.9) is a consequence of the following two inequalities:

(3.13)
(3.14)

which can be proved in the following way when .

If and (thus ), then it is easy to see that and

We shall prove so that . To this end, we consider the function

By considering the sign of the derivative , it is easy to see that the function achieves its minimum value at one of the two end points and , with

If , we then have . Consequently, for all and which implies

This proves (3.13). The inequality (3.14) can be proved in the same way. This completes the proof of (3.9) which further implies that and are analytic with respect to in the region enclosed by

The estimates (3.10)-(3.12) are simple consequences of Taylor’s theorem. ∎

To derive the representation of the numerical solution , we introduce some notations. Let be defined by

(3.15)
(3.16)
(3.17)

where we have set for ; this does not affect the value of , , upon solving (2.12). Similarly, we define

(3.18)
(3.19)
(3.20)

With these definitions, there are only a finite number of nonzero terms in the sequence , . Consequently, the generating function

is well defined (polynomial in ). Then, we have the following result.

Proposition 3.1.

For the time-stepping scheme (2.12), the semidiscrete solution can be represented by

(3.21)
(3.22)

where the operator is given by

(3.23)

with integration over the truncated contour defined in (3.7), oriented with increasing imaginary parts, with the parameters and satisfying the conditions of Lemma 3.3.

Proof.

In view of definition (2.10) and the identity (2.11), multiplying (2.12) by and summing up the results over yield

(3.24)

Then,

(3.25)

The function defined in (3.25) is analytic with respect to in a neighborhood of the origin. By Cauchy’s integral formula, it implies that for

where the second equality is obtained by the change of variables , with the contour defined in (3.8).

From Lemma 3.3, we see that both and are analytic with respect to in the region enclosed by , , and the two lines . Thus, is analytic with respect to . Because the values of on the two lines coincide, it follows that (by applying Cauchy’s integral formula)

(3.26)

where we have substituted (3.25) into the above equality and used the following (straightforward to check) identity in the last step:

with denoting the Laplace transform (in time) of the piecewise constant function .

Through the Laplace transform rule

(3.27)

one can derive (3.21) from (3.3). The proof of Proposition 3.1 is complete. ∎

3.4 Error estimate

In this subsection, we derive an error estimate for the numerical scheme (2.12). The following lemma is concerned with the difference between the kernels of (2.3) and (3.23). It will be used in the proof of Theorem 2.1.

Lemma 3.4.

Let be given and let be defined as in (2.9) with the parameters and satisfying the conditions of Lemma 3.3. Then, we have