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

 E∥ψ(⋅,tn)−ψn∥2L2(O)=O(τ1−αd/2)

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

 D(Δ)={ϕ∈H10(O):Δϕ∈L2(O)},

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

 ∂1−αtψ(x,t):=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩1Γ(α)∫t0(t−s)α−1∂ψ(x,s)∂sdsifα∈(0,1],1Γ(α−1)∫t0(t−s)α−2ψ(x,s)dsifα∈(1,2), (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

 ⎧⎪⎨⎪⎩∂tv(x,t)−Δ∂1−αtv(x,t)=f(x,t)(x,t)∈O×R+v(x,t)=0(x,t)∈∂O×R+v(x,0)=ψ0(x)x∈O (1.5)

plus the solution of the stochastic problem

 ⎧⎪⎨⎪⎩∂tu(x,t)−Δ∂1−αtu(x,t)=˙W(x,t)(x,t)∈O×R+u(x,t)=0(x,t)∈∂O×R+u(x,0)=0x∈O, (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

 E∥u(⋅,tn)−un∥2L2(O)=O(τ1−α/2−ϵ)

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

 E∥u(⋅,tn)−un∥2L2(O)=O(τ1−αd/2) α∈(0,2/d), \, d∈{1,2,3} (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])

 W(x,t)=∞∑j=1ϕj(x)Wj(t) (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])

 v(⋅,t)=∫t0E(t−s)f(⋅,s)ds, (2.2)

where the operator is given by

 E(t)ϕ:=12πi∫Γθ,κeztzα−1(zα−Δ)−1ϕdz∀ϕ∈L2(O), (2.3)

with integration over a contour on the complex plane,

 Γθ,κ ={z∈C:|z|=κ,|argz|≤θ}∪{z∈C:z=ρe±iθ,ρ≥κ} =:Γκθ,κ+Γθθ,κ. (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])

 u(⋅,t) =∫t0E(t−s)dW(⋅,s) (2.5) =∞∑j=1∫t0E(t−s)ϕjdWj(s). (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

 ¯∂1−ατun=1τ1−αn∑j=0bn−juj,n=0,1,2,…,N, (2.7)

where , , are the coefficients in the power series expansion

 (1−ζ)1−α=∞∑j=0bjζj. (2.8)

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

 δτ(ζ)=1−ζτfor ζ∈C∖[1,∞). (2.9)

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

 ˜v(ζ)=∞∑n=0vnζnfor ζ∈D (2.10)

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

 ˜v(eiθ)=limr→1−˜v(reiθ)

exists in . Then, we have

 ∞∑n=0(¯∂1−ατvn)ζn =∞∑n=01τ1−αn∑j=0bn−jvjζn (2.11) =(δτ(ζ))1−α∞∑j=0vjζj=(δτ(ζ))1−α˜v(ζ).

### 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

 un−un−1τ−Δ¯∂1−ατun=W(⋅,tn)−W(⋅,tn−1)τ. (2.12)

Equivalently, can be expressed as

 un= (1−ταb0Δ)−1un−1+ταn−1∑j=0bn−jΔ(1−ταb0Δ)−1uj (2.13) +(1−ταb0Δ)−1(W(⋅,tn)−W(⋅,tn−1)) = (1−ταb0Δ)−1un−1+ταn−1∑j=0bn−jΔ(1−ταb0Δ)−1uj +∞∑j=1(Wj(tn)−Wj(tn−1))(1+ταb0λj)−1ϕj.

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

 max1≤n≤NE∥u(⋅,tn)−un∥2L2(O)≤Cτ1−αd/2, (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 L2(Ω;L2(O))

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

 λj≥Cddd+2j2/d|O|−2/d (3.1)

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

Clearly, if for , then

 (1−ταb0Δ)−1un−1+ταn−1∑j=0bn−jΔ(1−ταb0Δ)−1uj∈L2(Ω;L2(O)). (3.2)

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

 (3.3)

In fact, we have

 E∥∥∥ℓ+m∑j=ℓ(Wj(tn)−Wj(tn−1))(1+ταb0λj)−1ϕj∥∥∥2 =Eℓ+m∑j=ℓ|Wj(tn)−Wj(tn−1)|2(1+ταb0λj)−2=ℓ+m∑j=ℓτ(1+ταb0λj)−2 ≤Cb−20τ1−2αℓ+m∑j=ℓj−4/d→0asℓ→∞.

Hence, for a fixed time step ,

 ℓ∑j=1(Wj(tn)−Wj(tn−1))(1+ταb0λj)−1ϕj,ℓ=1,2,…

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.
 ∞∑j=1(rαrα+λj)2≤Crαd/2∀r>0, (3.4) ∣∣∣1z+λj∣∣∣≤Cφ|z|+λj∀z∈Σφwithφ∈(0,π). (3.5)
###### Proof.

Clearly, Lemma 3.1 implies .

First, if ,

 ∞∑j=1(rαrα+λj)2≤∞∑j=1r2αCj4/d≤Cr2α≤Crαd/2,(d/2≤2 for d=1,2,3). (3.6)

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

 ∞∑j=1(rαrα+λj)2 ≤∞∑j=1(rαrα+Cj2/d)2 =M+1∑j=1(rαrα+Cj2/d)2+∞∑j=M+2(rαrα+Cj2/d)2=:I1+I2.

It is easy to see that and

 I2 ≤∫∞rαd/2(rαrα+Cs2/d)2ds=C−d/2rαd/2∫∞Cd/2(11+ξ2/d)2dξ≤Crαd/2,

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

 |z−ξ|sin(ω0)=|z|sin(ωξ)=λjsin(ωz).

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

 |z−ξ|≥|z|and|z−ξ|≥λj

which immediately implies

 |z−ξ|≥12(|z|+λj).

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

 |z−ξ|=|z|sin(ω0)sin(ωξ)≥|z|sin(φ)and|z−ξ|=λjsin(ω0)sin(ωz)≥λjsin(φ)

which immediately implies

 |z−ξ|≥sin(φ)2(|z|+λj).

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

 Γ(τ)θ,κ:={z∈Γθ,κ:|Im(z)|≤π/τ}. (3.7)

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

 Γ(τ)ρ:={z=−ln(ρ)/τ+iy:y∈R and |y|≤π/τ}. (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

 Γ(τ)ρ,\,\,\, Γ(τ)θ,κ,\,\,\,% and the two lines R±iπ/τwhenever0<κ≤min(1/T,−ln(ρ)/τ).

Furthermore, we have the following estimates:

 δτ(e−zτ)∈Σθ ∀z∈Γ(τ)θ,κ (3.9) C0|z|≤|δτ(e−τz)|≤C1|z| ∀z∈Γ(τ)θ,κ (3.10) |δτ(e−τz)−z|≤Cτ|z|2 ∀z∈Γ(τ)θ,κ (3.11) |δτ(e−τz)α−zα|≤Cτ|z|α+1 ∀z∈Γ(τ)θ,κ, (3.12)

where the constants , , and are independent of and .

###### Proof.

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

 0≤arg(1−e−zττ)≤arg(z) if0≤arg(z)≤θ, (3.13) if−θ≤arg(z)≤0, (3.14)

which can be proved in the following way when .

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

 cot(arg(1−e−τzτ)) =1−e−τ|z|cos(φ)cos(τ|z|sin(φ))e−τ|z|cos(φ)sin(τ|z|sin(φ)) =eτ|z|cos(φ)−cos(τ|z|sin(φ))sin(τ|z|sin(φ)) ≥1+τ|z|cos(φ)−cos(τ|z|sin(φ))sin(τ|z|sin(φ)) (Taylor's expansion) =1+ωcot(φ)−cos(ω)sin(ω)  (set ω=τ|z|sin(φ)∈[0,π)).

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

 g(ω):=1+ωcot(φ)−cos(ω)−sin(ω)cot(φ),ω∈[0,π].

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

 g(0)=0andg(π)=2+πcot(φ).

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

 0≤arg(1−e−τzτ)≤φ=arg(z).

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

 Γ(τ)ρ,\,\,\, Γ(τ)θ,κ\,\,\, % and the two lines R±iπ/τwhenever0<κ≤−ln(ρ)/τ.

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

 ¯∂τW(⋅,t0):=0 (3.15) ¯∂τW(⋅,t):=W(⋅,tn)−W(⋅,tn−1)τ fort∈(tn−1,tn],n=1,2,…,N (3.16) ¯∂τW(⋅,t):=0 for t>tN, (3.17)

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

 ¯∂τWj(t0):=0 (3.18) ¯∂τWj(t):=Wj(tn)−Wj(tn−1)τ fort∈(tn−1,tn],n=1,2,…,N (3.19) ¯∂τWj(t):=0 fort>tN. (3.20)

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

 ˜¯∂τW(⋅,ζ)=∞∑n=0¯∂τW(⋅,tn)ζn

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

 un =∫tn0Eτ(tn−s)¯∂τW(⋅,s)ds (3.21) =∞∑j=1∫tn0Eτ(tn−s)ϕj¯∂τWj(s)ds, (3.22)

where the operator is given by

 Eτ(t)ϕ:=12πi∫Γ(τ)θ,κeztzτezτ−1δτ(e−zτ)α−1(δτ(e−zτ)α−Δ)−1ϕdz∀ϕ∈L2(O) (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

 δτ(ζ)˜u(ζ)−δτ(ζ)1−αΔ˜u(ζ)=˜¯∂τW(⋅,ζ). (3.24)

Then,

 ˜u(ζ)=δτ(ζ)α−1(δτ(ζ)α−Δ)−1˜¯∂τW(⋅,ζ). (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

 un=12πi∫|ζ|=ρζ−n−1˜u(ζ)dζ=τ2πi∫Γ(τ)ρeztn˜u(e−zτ)dz,

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)

 un =τ2πi∫Γ(τ)ρeztn˜u(e−zτ)dz =τ2πi∫Γ(τ)θ,κeztn˜u(e−zτ)dz+τ2πi∫R+iπτeztn˜u(e−zτ)dz −τ2π%i∫R−iπτeztn˜u(e−zτ)dz =τ2πi∫Γ(τ)θ,κeztn˜u(e−zτ)dz =τ2πi∫Γ(τ)θ,κeztnδτ(e−zτ)α−1(δτ(e−zτ)α−Δ)−1˜¯∂τW(⋅,e−zτ)dz =τ2πi∫Γ(τ)θ,κeztnδτ(e−zτ)α−1(δτ(e−zτ)α−Δ)−1zezτ−1ˆ¯∂τW(⋅,z)dz, (3.26)

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

 ˜∂τW(⋅,e−zτ)=zezτ−1ˆ¯∂τW(⋅,z)

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

Through the Laplace transform rule

 L−1(ˆfˆg)(t)=∫t0L−1(ˆf)(t−s)L−1(ˆg)(s)ds, (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

 ∣∣∣zα−1(zα+λj)−1−zτ