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

TL;DR

This paper establishes sharp convergence rates for a time discretization scheme applied to stochastic time-fractional PDEs with white noise, demonstrating optimal error bounds and validating them through numerical experiments.

Contribution

The paper provides the first sharp-order error estimates for backward-Euler convolution quadrature applied to stochastic time-fractional PDEs with white noise, including optimal rates in one dimension.

Findings

Error estimate: $O( au^{1-rac{ ext{alpha} imes d}{2}})$ for the discretization

Optimal convergence rates for stochastic subdiffusion and diffusion-wave problems in 1D

Numerical examples confirm theoretical error bounds

Abstract

The stochastic time-fractional equation $\partial_t \psi -\Delta\partial_t^{1-\alpha} \psi = f + \dot W$ with space-time white noise $\dot W$ is discretized in time by a backward-Euler convolution quadrature for which the sharp-order error estimate \[ {\mathbb E}\|\psi(\cdot,t_n)-\psi_n\|_{L^2(\mathcal{O})}^2=O(\tau^{1-\alpha d/2}) \] is established for $\alpha\in(0,2/d)$, where $d$ denotes the spatial dimension, $\psi_n$ the approximate solution at the $n^{\rm th}$ time step, and $\mathbb{E}$ 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.

Tables

Table 1. Table 4.1. E ​ ( τ k ) 𝐸 subscript 𝜏 𝑘 E(\tau_{k}) and convergence rates in 1D.

$\alpha \backslash {\tau }_k$	$2^{-6}$	$2^{-7}$	$2^{-8}$	order
$\alpha =0.5$	1.075e-02	8.284e-03	6.382e-03	0.376 (0.375)
$\alpha =0.9$	2.825e-02	2.340e-02	1.921e-02	0.278 (0.275)
$\alpha =1.3$	6.340e-02	5.654e-02	5.004e-02	0.171 (0.175)
$\alpha =1.7$	1.415e-01	1.352e-01	1.275e-01	0.075 (0.075)

Table 2. Table 4.2. E ​ ( τ k ) 𝐸 subscript 𝜏 𝑘 E(\tau_{k}) and convergence rates in 2D.

$\alpha \backslash {\tau }_k$	$2^{-5}$	$2^{-6}$	$2^{-7}$	order
$\alpha =0.3$	3.848e-03	2.941e-03	2.335e-03	0.361 (0.35)
$\alpha =0.5$	9.992e-03	8.556e-03	7.395e-03	0.217 (0.25)
$\alpha =0.7$	2.115e-02	1.936e-02	1.762e-02	0.132 (0.15)
$\alpha =0.9$	3.997e-02	3.914e-02	3.851e-02	0.033 (0.05)

Equations

$$(\mathbb{E}\Vert \psi (\cdot ,t_n)-{\psi }_n{\Vert }_{L^2(\mathcal{O})}^2{)}^{\frac{1}{2}}=O({\tau }^{\frac{1}{2}-\frac{\alpha d}{4}})$$

$$\{\begin{array}{rlrl}& {\partial }_t\psi (x,t)-\Delta {\partial }_t^{1-\alpha }\psi (x,t)=f(x,t)+\dot{W}(x,t)& & (x,t)\in \mathcal{O}\times {\mathbb{R}}_{+}\\ & \psi (x,t)=0& & (x,t)\in \partial \mathcal{O}\times {\mathbb{R}}_{+}\\ & \psi (x,0)={\psi }_0(x)& & x\in \mathcal{O},\end{array}$$

$$D(\Delta )=\{\varphi \in H_0^1(\mathcal{O}):\Delta \varphi \in L^2(\mathcal{O})\},$$

$${\partial }_t^{1-\alpha }\psi (x,t):=\{\begin{array}{rlrl}& \frac{1}{\Gamma (\alpha )}{\int }_0^t(t-s{)}^{\alpha -1}\frac{\partial \psi (x,s)}{\partial s}\mathrm{d}s& & \text{mbox}if\alpha \in (0,1],\\ & \frac{1}{\Gamma (\alpha -1)}{\int }_0^t(t-s{)}^{\alpha -2}\psi (x,s)\mathrm{d}s& & \text{mbox}if\alpha \in (1,2),\end{array}$$

$${\partial }_t^{\alpha -1}{\partial }_t^{1-\alpha }\psi (x,t)=\{\begin{array}{rlrl}& \psi (x,t)-\psi (x,0)& & \text{mbox}if\alpha \in (0,1),\\ & \psi (x,t)& & \text{mbox}if\alpha \in (1,2),\end{array}$$

$${\partial }_t^{\alpha }\psi (x,t)-\Delta \psi (x,t)=\{\begin{array}{rlrl}& {\partial }_t^{\alpha -1}(f(x,t)+\dot{W}(x,t))-\Delta \psi (x,0)& & \text{mbox}if\alpha \in (0,1),\\ & f(x,t)+\dot{W}(x,t)-\Delta \psi (x,0)& & \text{mbox}if\alpha =1,\\ & {\partial }_t^{\alpha -1}(f(x,t)+\dot{W}(x,t))& & \text{mbox}if\alpha \in (1,2),\end{array}$$

$$\displaystyle\left\{\begin{array}[]{ll}\partial_{t}v(x,t)-\Delta\partial_{t}^{1-\alpha}v(x,t)=f(x,t)&(x,t)\in\mathcal{O}\times{\mathbb{R}}_{+}\\ v(x,t)=0&(x,t)\in\partial\mathcal{O}\times{\mathbb{R}}_{+}\\ v(x,0)=\psi_{0}(x)&x\in\mathcal{O}\end{array}\right.$$

$$\{\begin{array}{rlrl}& {\partial }_{t}u(x,t)-\Delta {\partial }_t^{1-\alpha }u(x,t)=\dot{W}(x,t)& & (x,t)\in \mathcal{O}\times {\mathbb{R}}_{+}\\ & u(x,t)=0& & (x,t)\in \partial \mathcal{O}\times {\mathbb{R}}_{+}\\ & u(x,0)=0& & x\in \mathcal{O}.\end{array}$$

$$(\mathbb{E}\Vert u(\cdot ,t_n)-u_n{\Vert }_{L^2(\mathcal{O})}^2{)}^{\frac{1}{2}}=O({\tau }^{\frac{1}{2}-\frac{\alpha }{4}-ϵ})$$

$$\displaystyle({\mathbb{E}}\|u(\cdot,t_{n})-u_{n}\|_{L^{2}(\mathcal{O})}^{2})^{\frac{1}{2}}=O(\tau^{\frac{1}{2}-\frac{\alpha d}{4}})\qquad\mbox{ $\alpha\in(0,2/d)$, \, $d\in\{1,2,3\}$}$$

$$W(x,t)=\sum _{j=1}^{\infty }{\varphi }_j(x)W_j(t)$$

$$v(\cdot ,t)={\int }_0^{t}E(t-s)f(\cdot ,s)\mathrm{d}s,$$

$$E(t)\varphi :=\frac{1}{2\pi \mathrm{i}}{\int }_{{\Gamma }_{\theta ,\kappa }}{e}^{zt}{z}^{\alpha -1}(z^{\alpha }-\Delta {)}^{-1}\varphi \mathrm{d}z\forall \varphi \in L^2(\mathcal{O}),$$

$${\Gamma }_{\theta ,\kappa }$$

$$=:{\Gamma }_{\theta ,\kappa }^{\kappa }+{\Gamma }_{\theta ,\kappa }^{\theta }.$$

$$u(\cdot ,t)$$

$$=\sum _{j=1}^{\infty }{\int }_0^{t}E(t-s){\varphi }_j\mathrm{d}{W}_j(s).$$

$${\bar{\partial }}_{\tau }^{1-\alpha }{u}_n=\frac{1}{{\tau }^{1-\alpha }}\sum _{j=0}^n{b}_{n-j}{u}_j,n=0,1,2,\dots ,N,$$

$$(1-\zeta {)}^{1-\alpha }=\sum _{j=0}^{\infty }{b}_j{\zeta }^j.$$

$${\delta }_{\tau }(\zeta )=\frac{1-\zeta }{\tau }\text{for }\zeta \in \mathbb{C}\backslash [1,\infty ).$$

$$\tilde{v}(\zeta )=\sum _{n=0}^{\infty }{v}_n{\zeta }^n\text{for }\zeta \in \mathbb{D}$$

$$\tilde{v}(e^{i\theta })=\underset{r\to 1_{-}}{\lim }\tilde{v}(r{e}^{i\theta })$$

$$\sum _{n=0}^{\infty }({\bar{\partial }}_{\tau }^{1-\alpha }{v}_n){\zeta }^n$$

$$=({\delta }_{\tau }(\zeta ){)}^{1-\alpha }\sum _{j=0}^{\infty }{v}_j{\zeta }^j=({\delta }_{\tau }(\zeta ){)}^{1-\alpha }\tilde{v}(\zeta ).$$

$$\frac{u_n-u_{n-1}}{\tau }-\Delta {\bar{\partial }}_{\tau }^{1-\alpha }{u}_n=\frac{W(\cdot ,t_n)-W(\cdot ,t_{n-1})}{\tau }.$$

$$u_n=$$

$$\displaystyle\qquad\qquad+\big{(}{\rm Id}-\tau^{\alpha}b_{0}\Delta\big{)}^{-1}\big{(}W(\cdot,t_{n})-W(\cdot,t_{n-1})\big{)}$$

$$=$$

$$\displaystyle\qquad\qquad+\sum_{j=1}^{\infty}\big{(}W_{j}(t_{n})-W_{j}(t_{n-1})\big{)}\big{(}1+\tau^{\alpha}b_{0}\lambda_{j}\big{)}^{-1}\phi_{j},$$

$$\displaystyle\max_{1\leq n\leq N}\,\bigg{(}{\mathbb{E}}\,\|u(\cdot,t_{n})-u_{n}\|_{L^{2}(\mathcal{O})}^{2}\bigg{)}^{\frac{1}{2}}\leq C\tau^{\frac{1}{2}-\frac{\alpha d}{4}},$$

$${\lambda }_j\ge \frac{C_{d}d}{d+2}{j}^{2/d}|\mathcal{O}{|}^{-2/d}$$

$$\displaystyle\big{(}{\rm Id}-\tau^{\alpha}b_{0}\Delta\big{)}^{-1}u_{n-1}+\tau^{\alpha}\sum_{j=0}^{n-1}b_{n-j}\Delta\big{(}{\rm Id}-\tau^{\alpha}b_{0}\Delta\big{)}^{-1}u_{j}\in L^{2}(\Omega;L^{2}(\mathcal{O})).$$

$$\displaystyle\sum_{j=1}^{\infty}\big{(}W_{j}(t_{n})-W_{j}(t_{n-1})\big{)}\big{(}1+\tau^{\alpha}b_{0}\lambda_{j}\big{)}^{-1}\phi_{j}\in L^{2}(\Omega;L^{2}(\mathcal{O})).$$

Full text

[TABLE]

is established for $\alpha\in(0,2/d)$, where $d$ denotes the spatial dimension, $\psi_{n}$ the approximate solution at the $n^{\rm th}$ time step, and $\mathbb{E}$ the expectation operator. In particular, the result indicates sharp 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.

1. Introduction

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

[TABLE]

where $\mathcal{O}\subset\mathbb{R}^{d}$, $d\in\{1,2,3\}$, denotes a bounded region with Lipschitz boundary $\partial\mathcal{O}$, $f(x,t)$ a given deterministic source function, $\psi_{0}(x)$ given deterministic initial condition, and $\dot{W}(x,t)$ a space-time white noise, i.e., the time derivative of a cylindrical Wiener process in $L^{2}(\mathcal{O})$. The underlying probability sample space for the stochastic noise is denoted by $\Omega$. The operator $\Delta:D(\Delta)\rightarrow L^{2}(\mathcal{O})$ denotes the Laplacian, defined on the domain

[TABLE]

and $\partial_{t}^{1-\alpha}\psi$ denotes the left-sided Caputo fractional time derivative of order $1-\alpha\in(-1,1)$, defined by (c.f. [18, pp. 91])

[TABLE]

where $\Gamma(s):=\int_{0}^{\infty}t^{s-1}e^{-t}{\rm d}t$ denotes Euler’s gamma function.

Problem (1.1) arises, e.g, when considering heat transfer in a material with thermal memory based on a modified Fick’s law [3, 10, 28, 32], subject to stochastic noise [4, 19, 30]. 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 $\partial_{t}^{\alpha}\psi(x,t):=\partial_{t}^{\alpha-1}\partial_{t}\psi(x,t)$ for $\alpha\in(1,2)$ and using the identity

[TABLE]

applying $\partial_{t}^{\alpha-1}$ to (1.1) yields another formulation of (1.1):

[TABLE]

where the case $\alpha=1$ 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

[TABLE]

plus the solution of the stochastic problem

[TABLE]

The stability and convergence of numerical solutions of (1.8) have been widely studied [2, 5, 26, 29, 31]. For example, if $f$ is smooth in time then numerical methods of up to order $6$ are available for approximating the solution of (1.8) and its equivalent formulations [14, 15, 16, 17, 21, 22, 26, 29]. In particular, the convolution quadrature generated by the backward Euler method yields a first-order convergence rate for solving (1.8).

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

[TABLE]

is proved [19, Remark 4.7, with $\rho=\alpha$]) for arbitrarily small $\epsilon>0$, where $u(\cdot,t_{n})$ and $u_{n}$ denote the mild solution and numerical solution of (1.9) at time $t_{n}$, respectively. The estimate (1.10) is “almost optimal” in the sense that the optimal approximation theoretic error estimate for functions having the regularity of the solution $u$ does not have the $\epsilon$ term in the exponent. We are not aware of any rigorous numerical analyses in the case $\alpha\in(0,1)$. In the case $\alpha=1$, error estimates for time discretization of the stochastic PDE (1.9) are proved in [12] and [1, 7, 11, 35] for Rothe’s method and the backward Euler method, respectively, with different spatial discretization methods. In particular, the convergence rate $O(\tau^{\frac{1}{4}})$ in time was proved, which corresponds to $\epsilon=0$ in (1.10). Some modern references on numerical analysis for stochastsic PDEs in the case $\alpha=1$ include [9, 13, 24, 36].

The error estimate (1.10) is consistent with the Hölder continuity of the solution $u\in C^{\gamma}([0,T];L^{2}(\Omega;L^{2}(\mathcal{O})))$ and the pathwise $\gamma$-Hölder continuity, with arbitrary $\gamma\in(0,\frac{1}{2}-\frac{\alpha d}{4})$; see Appendix A or [30, Corollary 1]. The aim of this paper is to prove, for general $d$-dimensional domains, the sharper order convergence estimate

[TABLE]

for time discretization of the stochastic PDE (1.9), where by “sharp” we mean that we are able to obtain an approximation theoretic convergence rate that is consistent with respect to the regularity of the solution in time. 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.9) 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 $C$, with/without a subscript, a generic constant independent of $n$ and $\tau$ 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 $\phi_{j}(x)$, $j=1,2,\dots$, denote the $L^{2}$-norm normalized eigenfunctions of the Laplacian operator $-\Delta$ corresponding to the eigenvalues $\lambda_{j}$, $j=1,2,\dots,$ arranged in nondecreasing order. The cylindrical Wiener process on $L^{2}(\mathcal{O})$ can be represented as (cf. [6, Proposition 4.7, with $Q=I$ and $U_{1}$ denoting some negative-order Sobolev space])

[TABLE]

with independent one-dimensional Wiener processes $W_{j}(t)$, $j=1,2,\dots$.

In the case $\psi_{0}=0$, the solution of the deterministic problem (1.8) can be expressed by (via Laplace transform, cf. [26, (3.11) and line 4 of page 12])

[TABLE]

where the operator $E(t):L^{2}(\mathcal{O})\rightarrow L^{2}(\mathcal{O})$ is given by

[TABLE]

with integration over a contour $\Gamma_{\theta,\kappa}$ on the complex plane,

[TABLE]

The angle $\theta$ above can be any angle such that $\pi/2<\theta<\min(\pi,\pi/\alpha)$ so that, for all $z$ to the right of $\Gamma_{\theta,\kappa}$ in the complex plane, $z^{\alpha}\in\Sigma_{\alpha\theta}:=\{z\in\mathbb{C}\backslash\{0\}:|\arg z|\leq\alpha\theta\}$ with $\alpha\theta<\pi$.

Correspondingly, the mild solution of (1.9) is define as (cf. [30] and [19, Proposition 2.7])

[TABLE]

This mild solution is well defined in $C^{\gamma}([0,T];L^{2}(\Omega;L^{2}(\mathcal{O})))$ for arbitrary $\gamma\in(0,\frac{1}{2}-\frac{\alpha d}{4})$; see Appendix A.

2.2. Convolution quadrature

Let $\{t_{n}=n\tau\}_{n=0}^{N}$ denote a uniform partition of the interval $[0,T]$ with a time step size $\tau=T/N$, and let $u^{n}=u(x,t_{n})$. Under the zero initial condition, the Caputo fractional time derivative $\partial_{t}^{1-\alpha}u(x,t_{n})$ can be discretized by the backward-Euler convolution quadrature [25] (also known as Grünwald-Letnikov approximation, cf. [5])

[TABLE]

where $b_{j}$, $j=0,1,2,\dots N$, are the coefficients in the power series expansion

[TABLE]

Here, $1-\zeta$ is the characteristic function of the backward-Euler method and we set

[TABLE]

For any given sequence $\{v_{n}\}_{n=0}^{\infty}\in\ell^{2}(L^{2}(\mathcal{O}))$, we denote

[TABLE]

which is referred to as the generating function of the sequence $\{v_{n}\}_{n=0}^{\infty}$ (see [26]). Clearly, $\widetilde{v}$ is an $L^{2}(\mathcal{O})$-valued analytic function in the unit disk ${\mathbb{D}}$ and the limit

[TABLE]

exists in $L^{2}(0,2\pi;L^{2}(\Omega))$. Then, we have

[TABLE]

2.3. Time-stepping scheme and main theorem

With the notations introduced in the last subsection, we discretize the fractional-order derivative $\partial_{t}^{1-\alpha}$ in (1.9) by using convolution quadrature in time to obtain

[TABLE]

Equivalently, $u_{n}$ can be expressed as

[TABLE]

where ${\rm Id}$ denotes the identity operator.

The main result of this paper is the following theorem.

Theorem 2.1**.**

Let $\alpha\in(0,2/d)$ with $d\in\{1,2,3\}$. Then, for each $n=1,2,\dots,N,$ the numerical solution $u_{n}$ given by (2.12) is well defined in $L^{2}(\Omega;L^{2}(\mathcal{O}))$ and converges to the mild solution $u(\cdot,t_{n})$ with sharp order of convergence, i.e., we have

[TABLE]

where ${\mathbb{E}}$ denotes the expectation operator and the constant $C$ is independent of $T$.

3. Proof of Theorem 2.1

3.1. The numerical solution in $L^{2}(\Omega;L^{2}(\mathcal{O}))$

In this subsection, we show that the numerical solution is well defined in $L^{2}(\Omega;L^{2}(\mathcal{O}))$. To this end, we use the following estimate for the eigenvalues of the Laplacian operator. For the simplicity of notations, we denote by $(\cdot,\cdot)$ and $\|\cdot\|$ the inner product and norm of $L^{2}(\mathcal{O})$, respectively.

Lemma 3.1** ([20, 23]).**

Let $\mathcal{O}$ denote a bounded domain in $\mathbb{R}^{d}$, $d\in\{1,2,3\}$. Suppose $\lambda_{j}$ denotes the $j^{\rm th}$ eigenvalue of the Dirichlet boundary problem for the Laplacian operator $-\Delta$ in $\mathcal{O}$. With $|\mathcal{O}|$ denoting the volume of $\mathcal{O}$, we have that

[TABLE]

for all $j\geq 1$, where $C_{d}=(2\pi)^{2}B_{d}^{-2/d}$ and $B_{d}$ denotes the volume of the unit $d$-dimensional ball.

Lemma 3.2**.**

Under the assumptions of Theorem 2.1, the numerical solution given by (2.13) is well defined in $L^{2}(\Omega;L^{2}(\mathcal{O}))$.

*Proof.  * Clearly, if $u_{j}\in L^{2}(\Omega;L^{2}(\mathcal{O}))$ for $j=0,\dots,n-1$, then

[TABLE]

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

[TABLE]

In fact, we have

[TABLE]

Hence, for a fixed time step $\tau$,

[TABLE]

is a Cauchy sequence in $L^{2}(\Omega;L^{2}(\mathcal{O}))$. Consequently, (3.3) is proved. In view of (2.13) and (3.2)-(3.3), the numerical solution $u_{n}$ is well defined in $L^{2}(\Omega;L^{2}(\mathcal{O}))$. ∎

3.2. A technical lemma

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

Lemma 3.3**.**

Let $\alpha\geq 0$ and $d\in\{1,2,3\}$. Then there exist constants $C$ and $C_{\varphi}$ such that

[TABLE]

where the constant $C$ depends on the dimension $d$ and the volume of the domain $\mathcal{O}$, and $C_{\varphi}$ depends on the angle $\varphi\in(0,\pi)$.

Proof.

Clearly, Lemma 3.1 implies $\lambda_{j}\geq Cj^{2/d}$, with a constant $C$ depending on the dimension $d$ and the volume of the domain $\mathcal{O}$.

First, if $0<r<1$,

[TABLE]

Second, if $r\geq 1$, by setting $M=\lfloor r^{\alpha d/2}\rfloor\geq 1$ to be the largest integer that does not exceed $r^{\alpha d/2}$, we have

[TABLE]

It is easy to see that $I_{1}\leq M+1\leq 2M\leq 2r^{\alpha d/2}$ and

[TABLE]

where the equality follows by changing the variable $s=C^{-d/2}r^{\alpha d/2}\xi$. This proves (3.4) in the case $r\geq 1$.

Finally, for the point $\xi=-\lambda_{j}+0\textrm{i}$ in the complex plane, we have $|\xi|=\lambda_{j}$ and $|z-\xi|=|z+\lambda_{j}|$. By looking at the triangle with three vertices $z$, [math], and $\xi$ with interior angles $\omega_{z}$, $\omega_{0}$, and $\omega_{\xi}$ at the three vertices, respectively, we have

[TABLE]

If $\omega_{0}\geq\pi/2$, then $|z-\xi|$ would be the length of the longest side of the triangle, i.e.,

[TABLE]

which immediately implies

[TABLE]

If $\omega_{0}\leq\pi/2$, then the angle condition $|{\rm arg}(z)|<\varphi$ implies $\omega_{0}>\pi-\varphi$. Hence, we have

[TABLE]

which immediately implies

[TABLE]

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

3.3. Solution representations

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

Let $\Gamma_{\theta,\kappa}^{(\tau)}$ denote the truncated piece of the contour $\Gamma_{\theta,\kappa}$ defined by

[TABLE]

For $\rho\in(0,1)$, let $\Gamma^{(\tau)}_{\rho}$ denote the segment of a vertical line defined by

[TABLE]

The following technical lemma will be used in this and next subsections, where ${\rm arccot}(\cdot)$ denotes the inverse of the cotangent function ${\rm cot}:(0,\pi)\rightarrow\mathbb{R}$.

Lemma 3.4**.**

Let $\alpha\geq 0$ and \theta\in\big{(}\frac{\pi}{2},{\rm arccot}(-\frac{2}{\pi})\big{)} be given, and let $\rho\in(0,1)$ be fixed, with $\delta_{\tau}(\zeta)$ defined in (2.9). Then, both $\delta(e^{-z\tau})$ and $(\delta(e^{-z\tau})^{\alpha}-\Delta)^{-1}$ are analytic with respect to $z$ in the region enclosed by

[TABLE]

Furthermore, we have the following estimates:

[TABLE]

where the constants $C_{0}$, $C_{1}$, and $C$ are independent of $\tau$ and $\kappa\in(0,\min(\frac{1}{T},-\frac{\ln(\rho)}{\tau}))$.

Proof.

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

[TABLE]

which can be proved in the following way when \frac{\pi}{2}\leq\theta\leq{\rm arccot}\big{(}-\frac{2}{\pi}\big{)}.

If ${\rm arg}(z)=\varphi\in[0,\theta]$ and $0\leq{\rm Im}(z)\leq\pi/\tau$ (thus $0\leq\tau|z|\sin(\varphi)\leq\pi$), then it is easy to see that {\rm arg}\big{(}\frac{1-e^{-\tau z}}{\tau}\big{)}\geq 0 and

[TABLE]

We shall prove \cot\big{(}{\rm arg}\big{(}\frac{1-e^{-\tau z}}{\tau}\big{)}\big{)}\geq\cot(\varphi) so that 0\leq{\rm arg}\big{(}\frac{1-e^{-\tau z}}{\tau}\big{)}\leq\varphi={\rm arg}(z). To this end, we consider the function

[TABLE]

whose derivative is

[TABLE]

In the case $\varphi\in(0,\frac{\pi}{2}]$, $g^{\prime}(\omega)\geq 0$ for $\omega\in[0,\pi]$. In the case $\varphi\in(\frac{\pi}{2},\pi)$, $g^{\prime}(\omega)\geq 0$ for $\omega\in[0,2\pi-2\varphi]$ and $g^{\prime}(\omega)\leq 0$ for $\omega\in[2\pi-2\varphi,\pi]$. In either case, the function $g(\omega)$ achieves its minimum value at one of the two end points $\omega=0$ and $\omega=\pi$, with

[TABLE]

If \frac{\pi}{2}\leq\theta\leq{\rm arccot}\big{(}-\frac{2}{\pi}\big{)}, we then have $g(\pi)\geq 0$. Consequently, $g(\omega)\geq 0$ for all $\omega\in[0,\pi]$ and \cot\big{(}{\rm arg}\big{(}\frac{1-e^{-\tau z}}{\tau}\big{)}\big{)}\geq\cot(\varphi) which implies

[TABLE]

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 $\delta(e^{-z\tau})$ and $(\delta(e^{-z\tau})^{\alpha}-\Delta)^{-1}$ are analytic with respect to $z$ in the region enclosed by

[TABLE]

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

Remark 3.1**.**

The condition $\kappa\leq\frac{1}{T}$ is not needed in the proof of this lemma, but is needed in the estimates of the next subsection such as (3.4).

To derive the representation of the numerical solution $u_{n}$, we introduce some notations. Let $\bar{\partial}_{\tau}W$ be defined by

[TABLE]

where we have set $\bar{\partial}_{\tau}W(\cdot,t)=0$ for $t>t_{N}$; this does not affect the value of $u_{n}$, $n=1,2,\dots,N$, upon solving (2.12). Similarly, we define

[TABLE]

With these definitions, there are only a finite number of nonzero terms in the sequence $\bar{\partial}_{\tau}W(\cdot,t_{n})$, $n=0,1,2,\dots$. Consequently, the generating function

[TABLE]

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

Proposition 3.2**.**

For the time-stepping scheme (2.12), the semidiscrete solution $u_{n}$ can be represented by

[TABLE]

where the operator $E_{\tau}(\cdot)$ is given by

[TABLE]

with integration over the truncated contour $\Gamma_{\theta,\kappa}^{(\tau)}$ defined in (3.7), oriented with increasing imaginary parts, with the parameters $\kappa$ and $\theta$ satisfying the conditions of Lemma 3.4.

Proof.

In view of definition (2.10) and the identity (2.11), multiplying (2.12) by $\zeta^{n}$ and summing up the results over $n=0,1,2,\dots$ yield

[TABLE]

Then,

[TABLE]

The function $\widetilde{u}(\zeta)$ defined in (3.25) is analytic with respect to $\zeta$ in a neighborhood of the origin. By Cauchy’s integral formula, it implies that for $\rho\in(0,1)$

[TABLE]

where the second equality is obtained by the change of variables $\zeta=e^{-z\tau}$, with the contour $\Gamma^{(\tau)}_{\rho}$ defined in (3.8).

From Lemma 3.4, we see that both $\delta(e^{-z\tau})$ and $(\delta(e^{-z\tau})^{\alpha}-\Delta)^{-1}$ are analytic with respect to $z$ in the region $\Sigma\subset\mathbb{C}$ enclosed by $\Gamma^{(\tau)}_{\rho}$, $\Gamma^{(\tau)}_{\theta,\kappa}$, and the two lines $\mathbb{R}\pm\textrm{i}\pi/\tau$. Thus, $e^{zt_{n}}\widetilde{u}(e^{-z\tau})$ is analytic with respect to $z\in\Sigma$. Because the values of $e^{zt_{n}}\widetilde{u}(e^{-z\tau})$ on the two lines $\mathbb{R}\pm\textrm{i}\pi/\tau$ coincide, it follows that (by applying Cauchy’s integral formula)

[TABLE]

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

[TABLE]

with $\widehat{\bar{\partial}_{\tau}W}$ denoting the Laplace transform (in time) of the piecewise constant function $\bar{\partial}_{\tau}W$.

Through the Laplace transform rule

[TABLE]

one can derive (3.21) from (3.3). The proof of Proposition 3.2 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.5**.**

Let $\alpha\in(0,2/d)$ be given and let $\delta_{\tau}(\zeta)$ be defined as in (2.9) with the parameters $\kappa$ and $\theta$ satisfying the conditions of Lemma 3.4. Then, we have

[TABLE]

Proof.

By the triangle inequality and Lemma 3.4, we have

[TABLE]

where

[TABLE]

In the estimates above we have used the following inequality (cf. [8, inequality C.1])

[TABLE]

The last inequality is due to Lemma 3.3 together with the angle condition ${\rm arg}(z^{\alpha})\leq\alpha\theta<\pi$ and ${\rm arg}(\delta_{\tau}(e^{-\tau z})^{\alpha})\leq\alpha\theta<\pi$ (cf. Lemma 3.4). Furthermore, we have

[TABLE]

The proof of Lemma 3.5 is complete. ∎

Now, we start to prove Theorem 2.1. From (2.6) and (2.3) we see that the mild solution admits the decomposition

[TABLE]

with

[TABLE]

Also, (3.21) and (3.23) imply

[TABLE]

with

[TABLE]

Comparing (3.29) and (3.32) yields

[TABLE]

Then Theorem 2.1 is a consequence of the following lemma. The proof of Theorem 2.1 is complete. ∎

Lemma 3.6**.**

Under the assumptions of Theorem 2.1, we have ${\mathcal{E}}_{\tau}(t_{n}),{\mathcal{G}}_{\tau}(t_{n}),{\mathcal{H}}_{\tau}(t_{n})\in L^{2}(\Omega;L^{2}(\mathcal{O}))$, satisfying the following estimate:

[TABLE]

*Proof.  * First, we estimate ${\mathcal{H}}_{\tau}(t_{n})$. By choosing a number $\beta\in(\alpha d/2,1)$ and using Lemma 3.3, we have

[TABLE]

Next, we estimate ${\mathcal{E}}_{\tau}(t_{n})$. To this end, we apply Lemma 3.5 and obtain

[TABLE]

By using the expression of $\mathcal{E}_{\tau}(t_{n})$, we have

[TABLE]

where the last inequality follows from Lemma 3.3. Since $t_{n}\geq\tau$ and

[TABLE]

by choosing $\kappa\leq\frac{2}{t_{n}|\sin(\theta)|}$, we have

[TABLE]

Finally, we estimate $\mathcal{G}_{\tau}(t_{n})$. Because $\bar{\partial}_{\tau}W_{j}(t_{n})=\frac{1}{\tau}\int_{t_{n-1}}^{t_{n}}{\rm d}W_{j}(s)$, we obtain

[TABLE]

Then,

[TABLE]

By using the expression (3.33), for $|s-\xi|\leq\tau$ we have

[TABLE]

where we have used $|1-e^{z(s-\xi)}|\leq C|z(s-\xi)|\leq C|z|\tau$ and (3.28) in deriving the second to last inequality. Substituting the last inequality into (3.4) yields

[TABLE]

where the last inequality can be estimated in the same way as (3.4). ∎

4. Numerical examples

In this section, we present three numerical examples to illustrate the theoretical analyses.

Example 1. We first consider the one-dimensional stochastic time-fractional equation

[TABLE]

for $0\leq x\leq 1$, $0<t\leq 1$, with homogenous Dirichlet boundary condition and [math] initial condition. In the above equation,

[TABLE]

$\varepsilon$ is a given constant, and $W$ the cylindrical Wiener process. In the absence of white noise, the exact solution would be $u_{d}(x,t)=t^{2}x^{2}(1-x)^{2}$, which corresponds to the exact mean of the stochastic solution.

We discretize the problem (4.1) in time by using the scheme (2.12) and, in space, by continuous piecewise linear finite element method. Here, $h=1/M$ denotes the spatial mesh size and $U^{n}(x)$ the numerical solution of the fully discrete scheme. We take $\tau=h=2^{-5}$ and $\varepsilon=0.1$. For each computation, $I=1000$ independent realizations are performed with different Wiener processes. For each realization $\omega_{i}$, $i=1,\dots,I$, we generate $M$ independent Brownian motions $W_{j}(t)$, $j=1,\dots,M$. In Figure 4.1(left), we present the exact solution $u_{d}$ of the deterministic problem, the mean value of numerical solutions for (4.1), and the standard deviation, respectively, at $t_{n}=1$. Moreover, the numerical approximations $U^{n}(x,\omega_{i})$, $i=1,2,3$ of $u(x,t_{n},\omega_{i})$, with three independent realizations, are given in Figure 4.1(right) at $t_{n}=1$. The numerical simulations in Figure 4.1 are performed by taking $\alpha=0.5$. Similar results are shown in Figure 4.2 for $\alpha=1.3$. Because the solution has $C^{\min(\frac{1}{\alpha}-\frac{1}{2},1)}(\overline{\Omega})$ pathwise regularity in space (cf. [30, Proposition 2]), the numerical solution for $\alpha=0.5$ (the solution is $C^{1}(\overline{\Omega})$) is smoother than the numerical solution for $\alpha=1.3$ (the solution is $C^{0.27}(\overline{\Omega})$); see Figures 4.1 and 4.2 for a visual comparison.

Example 2. We next consider the convergence rate of the numerical scheme (2.12) for (4.1) with $\varepsilon=1$. The problem (4.1) is discretized using backward-Euler convolution quadrature and a linear Galerkin finite element method. To investigate the convergence rate, we consider $I=1000$ independent realizations for each time step $\tau_{k}=2^{-k}$, $k=5,\dots,8$. In order to focus on the time discretization error, we solve the time-discrete stochastic PDE (2.12) using a sufficiently small spatial mesh size $h=1/M=2^{-9}$ so that the spatial discretization error is relatively negligible. Then the error $E(\tau_{k})$ is computed by

[TABLE]

for $k=6,7,8$.

In [26], it is proved that the backward-Euler convolution quadrature for time-fractional PDE (1.8) is first-order convergent. Thus, by Theorem 2.1, the convergence order of the scheme (2.12) for problem (4.1) should be $O(\tau^{\frac{1}{2}-\frac{\alpha}{4}})$ in a one-dimensional spatial domain. Consequently, we expect the error $E(\tau_{k})$ to have the convergence rate

[TABLE]

for successive halvings of the time step. We test the above result by taking $\alpha=0.5$ and $0.9$ for a subdiffusion setting and $\alpha=1.3$ and $1.7$ for a diffusion-wave setting. From (4.3), $\log_{2}\frac{E(\tau_{k})}{E(\tau_{k+1})}\approx 0.375$ for $\alpha=0.5$, $\log_{2}\frac{E(\tau_{k})}{E(\tau_{k+1})}\approx 0.275$ for $\alpha=0.9$, $\log_{2}\frac{E(\tau_{k})}{E(\tau_{k+1})}\approx 0.175$ for $\alpha=1.3$, and $\log_{2}\frac{E(\tau_{k})}{E(\tau_{k+1})}\approx 0.075$ for $\alpha=1.7$. Clearly, the results in Table 4.1 ($t_{n}=1$) illustrate the sharp convergence rate.

Example 3. Lastly, we consider the stochastic time-fractional equation

[TABLE]

in the two-dimensional spatial domain $[0,1]\times[0,1]$, with homogenous Dirichlet boundary condition and [math] initial condition. Here, we choose $0<t\leq 1$ and

[TABLE]

for $x=(x_{1},x_{2})\in[0,1]\times[0,1]$. The exact solution of the corresponding deterministic problem is $u_{d}(x,t)=t^{2}x_{1}^{2}x_{2}^{2}(1-x_{1})^{2}(1-x_{2})^{2}$.

We solve the stochastic equation (4.4) using the backward-Euler scheme (2.12), where spatial discretization is effected by the standard piecewise linear Galerkin finite element method. A uniform triangular partition with 50 nodes in each direction is used. Similarly, we choose $\tau_{k}=2^{-k}$, $k=4,\dots,7$, and consider $I=1000$ independent realizations to investigate the temporal convergence rate. For each realization $\omega_{i}$, $i=1,\dots,I$, we generate $2500$ independent Brownian motions $W_{j}(t)$, $j=1,\dots,2500$. The mesh size $h=\frac{\sqrt{2}}{50}$ is fixed so that spatial error is relatively negligible. Then, the error (4.2) is computed for each fixed time step $\tau_{k}$ and presented in Table 4.2 at $t_{n}=1$. By Theorem 2.1, the convergence rate of the scheme (2.12) for problem (4.4) is $O(\tau^{\frac{1}{2}-\frac{\alpha}{2}})$ in the two spatial dimensional setting. Clearly, the numerical results are consistent with the theoretical analyses given in Theorem 2.1.

5. Conclusion

We considered the stability and convergence of numerical approximations of a stochastic time-fractional PDE by using the backward Euler convolution quadrature in time. By means of a discrete analogue of the inverse Laplace transform, we derived an integral representation of the numerical solution which was then used to prove the sharp convergence rate of the numerical approximation.

Instead of the contour

[TABLE]

used in [26], we have used the contour $\Gamma_{\theta,\kappa}^{(\tau)}$ given in (3.7) for the analysis in this paper. The contour $\Gamma_{\theta,\kappa}^{(\tau)}$ excludes the origin and thus can be used not only for the Dirichlet Laplacian but also for the Neumann Laplacian (whose spectrum includes the origin). Similarly, for the Dirichlet Laplacian, it is not necessary to show that $\delta_{\tau}(e^{-z\tau})\in\Sigma_{\theta}$ as in (3.9). Using this result, the analyses in this paper can be naturally extended to the Neumann Laplacian.

The paper focuses on semidiscretization in time by convolution quadrature. The main contribution of the paper is to show the possibility of removing the $\epsilon$ term in the previously obtained error estimate (1.10) and to establish a foundation for the further analysis of spatial discretization considered in [8]. Our analysis is based on specific growth properties of the eigenvalues of the Laplacian operator (cf. Lemma 3.1), and thus may not be directly extended to more general abstract operators such as semilinear problems and multiplicative space-time white noises. For example, for the multiplicative noise case, the identities (3.29)-(3.33) do not hold and thus the analysis becomes more complicated. Thus, extension to semilinear problems and multiplicative space-time white noises remains open and certainly should be a subject of future research.

Theorem 2.1 can be extended to higher order moments by using the Burkhölder-Davis-Gundy inequality (cf. [33, (1.1)], [6, (6.29)], or [27])$:$ for all $p\in(1,\infty)$ there exists $C_{p}>0$ such that

[TABLE]

where $L_{2}^{0}$ denotes the space of Hilbert-Schmidt operators on $L^{2}(\mathcal{O})$. Let

[TABLE]

which is a Hilbert-Schmidt operator satisfying (see (3.31))

[TABLE]

By using the definition of $\mathcal{H}_{\tau}(t_{n})$ in (3.4), we have

[TABLE]

where the last inequality utilizes the result of (3.4). Similarly, the following estimates can be proved by using (3.4) and (3.4):

[TABLE]

Substituting these estimates into (3.4) yields the $p$-moment estimate:

[TABLE]

The estimate above also implies the existence of a random variable $\mathcal{C}$ having finite moments of any order independent of $\tau$ such that the following pathwise estimate holds:

[TABLE]

The convergence rate proved in this article is optimal with respect to the regularity of the solution in time. However, whether it is the highest possible convergence rate among all possible numerical methods is unknown. For stochastic ODEs, the convergence rates of some numerical methods (such as Milstein’s method) may be higher than the regularity of the solution in time. But Milstein’s method does not yield higher convergence rate for the stochastic PDE problem considered in this article. For example, in the case $\alpha=d=1$, Milstein’s method (equivalent to Euler Maruyama with additive noise)

[TABLE]

does not converge at all (as it is an explicit scheme, which requires a CFL condition that cannot be satisfied by a semi-discretization in time). Even if we modify the Milstein’s method to be an implicit scheme

[TABLE]

the scheme only has strong convergence rate $O(\tau^{1/4})$. This is different from stochastic ODEs, for which the strong convergence rates of Milstein’s method and the corresponding implicit scheme are $O(\tau)$ (higher than the temporal regularity of the solution). This difference between stochastic PDEs and stochastic ODEs is due to the fact that $\Delta u$ does not have the same temporal regularity as $u$ when the PDE is driven by a space-time white noise (2.1). In particular, $\Delta u$ is not Hölder continuous in time in the space $L^{2}(\mathcal{O})$ (but Milstein’s method requires Hölder continuity of $\Delta u$ to achieve a higher convergence rate). In the case of colored noise (instead of space-time white noise), the implicit Euler scheme (5.6) can achieve a better convergence rate than the temporal Hölder regularity of the solution; see [34].

Appendix A The mild solution in $C^{\gamma}([0,T];L^{2}(\Omega;L^{2}(\mathcal{O})))$

In the cases $\alpha<\min(1,2/d)$ and $\alpha>1$ with $d=1$, the mild solution of (1.9) (with space-time white noise) has been studied in different function spaces under different settings. For example, see [19, 30]. For the reader’s convenience, in this appendix, we illustrate that the mild solution given by (2.5)-(2.6) is indeed well defined in $C^{\gamma}([0,T];L^{2}(\Omega;L^{2}(\mathcal{O})))$, a result used for the numerical analysis in this paper.

Theorem A.1**.**

The mild solution defined by (2.5)-(2.6) is in $C^{\gamma}([0,T];L^{2}(\Omega;L^{2}(\mathcal{O})))$ for arbitrary $\gamma\in(0,\frac{1}{2}-\frac{\alpha d}{4})$.

Proof.   In (2.6), the formula (2.3) implies

[TABLE]

for a deterministic time-independent function $\phi_{j}\in L^{2}(\mathcal{O})$ and with the deterministic space-independent function $h_{j}(\cdot)$ given by

[TABLE]

By the theory of the Ito integral and the identity (A.1), each term in (2.6) is well defined in $C([0,T];L^{2}(\Omega;L^{2}(\mathcal{O})))$. Because the one-dimensional Wiener processes $W_{j}(s)$, $j=1,2,\dots$, are independent of each other, it follows that

[TABLE]

where $\beta\in(1,2-\alpha d/2)$.

In view of Lemma 3.3, \displaystyle\sum_{j=1}^{\infty}\bigg{|}\frac{r^{\alpha}}{r^{\alpha}+\lambda_{j}}\bigg{|}^{2}\leq Cr^{\alpha d/2} implies \displaystyle\sum_{j=\ell}^{\infty}\bigg{|}\frac{r^{\alpha}}{r^{\alpha}+\lambda_{j}}\bigg{|}^{2}\rightarrow 0 as $\ell\rightarrow\infty$ and

[TABLE]

The Lebesgue dominated convergence theorem implies that

[TABLE]

Similarly, \displaystyle\sum_{j=\ell}^{\infty}\bigg{|}\frac{\kappa^{\alpha}}{\kappa^{\alpha}+\lambda_{j}}\bigg{|}^{2}\leq C\kappa^{\alpha d/2} implies \displaystyle\sum_{j=\ell}^{\infty}\bigg{|}\frac{\kappa^{\alpha}}{\kappa^{\alpha}+\lambda_{j}}\bigg{|}^{2}\rightarrow 0 as $\ell\rightarrow\infty,$ and

[TABLE]

Again, the Lebesgue dominated convergence theorem implies that

[TABLE]

Overall, we have

[TABLE]

which implies that the sequence

[TABLE]

is a Cauchy sequence in $C([0,T];L^{2}(\Omega;L^{2}(\mathcal{O})))$. Consequently, the sequence converges to a function $u\in C([0,T];L^{2}(\Omega;L^{2}(\mathcal{O})))$, which is the mild solution defined by (2.6).

Let $L_{2}^{0}$ denote the space of Hilbert-Schmidt operators on $L^{2}(\mathcal{O})$ (cf. [6, Appendix C]) with the operator norm

[TABLE]

The above analysis clearly shows that

[TABLE]

In view of [6, Proposition 4.20 and page 99], the stochastic integral (2.5) is well defined in $L^{2}(\Omega;L^{2}(\mathcal{O}))$, and (2.5) coincides with the series representation (2.6) ([6, section 4.2.2]).

Similar to the estimate (A), by considering

[TABLE]

we have

[TABLE]

Because \big{|}\frac{e^{zh}-1}{h^{\gamma}}\big{|}\leq C|z|^{\gamma} on the contour $\Gamma_{\theta,\kappa}$ (on which ${\rm Re}(z)\leq 0$ when $|z|\geq\kappa$), it follows that

[TABLE]

where the second last inequality requires $\beta>1$ for the improper integral $\int_{\Gamma_{\theta,\kappa}^{\kappa}}\frac{1}{|z|^{\beta}}|{\rm d}z|$ to be convergent (then the estimates follow similarly as for (A)-(A)), and the last inequality requires $\alpha d/2+\beta+2\gamma-3<-1$. This requires $2\gamma<1-\frac{\alpha d}{2}$. Also,

[TABLE]

which again requires $\beta>1$ and $2\gamma<1-\frac{\alpha d}{2}$ for the convergence of the improper integrals.

Substituting the last two results into (A) yields that $u\in C^{\gamma}([0,T];L^{2}(\Omega;L^{2}(\mathcal{O})))$ for arbitrary $\gamma\in(0,\frac{1}{2}-\frac{\alpha d}{4})$. ∎

Acknowledgements. The authors would like to thank Professor Xiaobing Feng for many valuable suggestions and comments.

Funding. The research of M. Gunzburger and J. Wang was supported in part by the USA National Science Foundation grant DMS-1315259 and by the USA Air Force Office of Scientific Research grant FA9550-15-1-0001. The work of B. Li was supported in part by the Hong Kong RGC grant 15300817.