An efficient numerical approach for stochastic evolution PDEs driven by random diffusion coefficients and multiplicative noise

1.   Introduction

Stochastic partial differential equations (SPDEs) appears in many fields of science and engineering, and have been subject of many theoretical and numerical investigations. It is commonly believed that incorporating noise and/or uncertainty into models is closer to reality in mathematical modeling, due to the existence of uncertainty stemming from various sources such as thermal fluctuation, impurities of materials and so on. As an active area of research, numerical study of stochastic evolution equations (SEEs) has attracted increasing attention in the past decades; see, e.g., monographs ^{[33,40,45,53,56,70]} and references therein. Although much progress has been made, it is still far from being satisfactory due to the numerical approximations of SEEs, especially for the SEEs with non-globally Lipschitz nonlinearity, encounter all the difficulties that may arise in solving deterministic differential equations on one hand, and caused by the effect of nonlinearity, infinite dimensional operator and driving noise on the other hand, see, e.g., ^{[7,11,13,14,15,16]} and references therein. The present work focus on the SEEs perturbed by a smooth random diffusion coefficient field as well as multiplicative force noise, and aims to propose and analyze an efficient numerical method for this equation.

When considering the numerical approaches for SEEs with various noises, two categories of convergence errors may be involved, namely weak error and strong error. The former is related to the approximation of the probability law of the solution. Concerning weak convergence error of numerical methods for SEEs, we refer to, for instance, ^{[2,6,8,12,16,19,20,21,27,32,41,42,50,58,66]} and references therein for a list of literature in this direction. Unlike weak convergence error, the strong convergence error measures the deviation from the trajectory of an exact solution. It has been extensively investigated in various types of SPDEs, see, e.g., ^{[1,5,7,9,11,13,14,15,18,24,26,29,30,34,35,36,39,43,44,51,52,59,63,64,65,69]} and references therein. We mention here some works on strong convergence of the numerical schemes for linear SEEs with additive or multiplicative noise. For example, Allen et al. ^[1] described, analyzed and compared the finite element and difference methods for parabolic SPDEs driven by additive white noise. Du et al. ^[24] investigated numerical solutions of linear SEEs perturbed by special additive noises, ranging from the space time white noise to colored noises generated by some infinite dimensional Brownian motions with a prescribed covariance operator. Yan ^[69] studied the finite element method for linear SEEs with multiplicative noise in multidimensional case. The case of strong convergence of nonlinear SEEs is generally more subtle and challenging, and has received widely attention in the research community in recent years. For instance, Kloeden et al. ^[35,39] proposed a discretization based on the Galerkin method in space and exponential integrator in time for the nonlinear SEEs with cylindrical additive noise. Kruse ^[44] analyzed the strong convergence error for a finite element method/linear implicit Euler spatio-temporal discretization of semilinear SEEs with multiplicative noise and Lipschitz continuous nonlinearities, and deduced the optimal error estimates. Wang ^[64] derived strong convergence results for a spatio-temporal discretization of the semilinear SEEs with additive noise, where the approximation in space was performed by a standard finite element method and in time by a linear implicit Euler method. Moreover it was shown how exactly the strong convergence rate of the full discretization relies on the regularity of the driven process. Kov$\acute{\mathrm{a}}$cs et al. ^[43] used Euler type splitstep method to study the semidiscretisation in time of the stochastic Allen-Cahn equation perturbed by smooth additive Gaussian noise, and showed that the strong convergence rate is 1/2 with respect to the step size. Cui et al. ^[11,13,14] analyzed and obtained the strong convergence rate of the finite difference approximation and splitting scheme of the stochastic nonlinear Schrödinger equation. Bréhier et al. ^[7] derived the optimal strong convergence rate for the explicit splitting schemes of the stochastic Allen-Cahn equation. Liu et al. ^[52] proposed a general theory of optimal strong error estimation for some drift-implicit Euler schemes of a second-order nonlinear SPDE with monotone drift driven by a multiplicative infinite-dimensional Wiener process.

In this paper, we consider the SEEs with both multiplicative force noise and random diffusion coefficient field, which has not yet been addressed in the literature to the best of our knowledge. The main contributions/novelties of this paper are as follows:

$\bullet$ The well-posedness of the considered stochastic equation is established. That is, the existence, uniqueness, and stability of the mild solution is proved.

$\bullet$ The diffusion coefficient considered in the current work is a bounded log-Whittle-Matérn Gaussian random field with a parametrized covariance function whose regularity can be controlled by a parameter. Therefore different cases can be tested and compared in a convenient way.

$\bullet$ A sampling approach called approximation circulant embedding with padding ^[23,57,67] is employed to render the equation solvable. Then for each sample diffusion coefficient, a time-stepping scheme based on a semi-implicit Euler-Maruyama approach is constructed for the resulting equation. The standard piecewise linear finite element method is employed for the spatial discretization. The main theoretical result is the proof of the strong convergence rate $\mathcal{O}(h^{2-\varepsilon_0}+\Delta t^{\frac{1}{2}})$ of the full discretization under certain assumptions, where $\varepsilon_0$ is an infinitesimal positive number, $h$ and $\Delta t$ are respectively the spatial and temporal mesh sizes.

The paper is organised as follows. In Section 2, we establish the well-posedness of the considered problem under given assumptions. The sampling method for the random diffusion coefficient field as well as the spatio-temporal full discretization are presented in Section 3. We devote to deriving the strong error estimate of the proposed fully discrete scheme by using semigroup approach and the stochastic calculus tools in Section 4, and validate the theoretical results by numerical experiments in Section 5.

2.   Problem and its well-posedness

We start by defining our problem. Let $T > 0$, $D\in\mathbb{R}^d$, $d\in\{1, 2, 3\}$, be a bounded open spatial domain with smooth boundary. To be specific, we consider $D: = (0, 1)^d$ in this work. Let $L^2(D)$ and $H_0^{\gamma}(D)$ be classical Sobolev spaces, $\gamma\ge0$. $\mathcal{L}(L^2(D))$ represents the space of bounded linear operators $A$: $L^2(D)\to L^2(D)$ equipped with operator norm $\rVert A\rVert_{_{{\mathcal L}(L^2(D))}} = \sup\limits_{u\ne0}\frac{\rVert Au\rVert_{_{L^2(D)}}}{\rVert u\rVert_{_{L^2(D)}}}$. $(\Omega, \mathcal F, \{\mathcal F_t\}_{t\ge 0}, \mathbb{P})$ is a filtered probability space with a normal filtration $\{\mathcal{F}_t\}_{t\ge0}$. Additionally, we denote $v(x, \omega)\in L^2(\Omega, L^2(D))$ if

where the norm $\rVert v \rVert_{_{L^2(\Omega, L^2(D))}}$ is defined by

with $\mathbb{E}[\cdot]$ being the expectation in the probability space $(\Omega, \mathcal{F}, \mathbb{P})$. $L^2(\Omega, L^2(D))$ is also known as the space of the mean-square integrable random variables. Let $W(t, x)$ be a $\mathcal{F}_t$-adapted $H_0^{\gamma}(D)$-valued Wiener process with covariance operator $Q$, where $Q$ is a positive definite and symmetric operator with orthonormal eigenfunctions {$\phi_j(x)\in H_0^{\gamma}(D):j\in\mathbb{N}$} and corresponding positive eigenvalues $\{q_j\}$; see, e.g., ^[44,64,69] for more details.

Let $Q^{\frac{1}{2}}(H_0^{\gamma}(D)): = \{Q^{\frac{1}{2}}v: v\in H_0^{\gamma}(D)\}$. Let ${\mathcal L}_Q$ be the set of linear operators $B: Q^{\frac{1}{2}}(H_0^{\gamma}(D))\to L^2(D)$, which satisfies

${\mathcal L}_Q$ endowed with the norm $\rVert B \rVert_{{\mathcal L}_Q}: = \Big(\sum_{j = 1}^\infty \rVert BQ^\frac{1}{2}\phi_j \rVert_{_{L^2(D)}}^2\Big)^\frac{1}{2}$ is the space of Hilbert-Schmidt operators ^[28]. We will also use the space $L^2(\Omega, {\mathcal L}_Q)$ of all random Hilbert-Schmidt operators $B: \Omega\to {\mathcal L}_Q$, equipped with the norm

Throughout the paper we use $c$, with or without subscripts, to mean generic positive constants (independent of $\omega$ in particular), which may not be the same at different occurrences.

Our point of interest is the SEE with random diffusion coefficient and multiplicative noise, written in the abstract form:

where $L: = -\nabla\cdot(a(x, \omega)\nabla)$ is the elliptic operator with the coefficient $a(x, \omega)$ being a bounded log-Gaussian random field with the scale parameter $\mathscr{ \varepsilon}$, i.e., there exists two constants $a_{\min}$ and $a_{\max}$ such that: for almost every $x\in\bar{D}$ and $\omega\in\Omega$,

Clearly the satisfaction of (2.3) relies on the uniform boundedness of the Gaussian random variable $z(x, \omega)$. Notice that (2.3) does not hold if $z(x, \omega)$ is Gaussian random variable without any restriction. Here we assume that for any $x\in\bar{D}$, $z(x, \omega)$ is a truncated Gaussian random variable ^[10,37], such that $z_{\min}\leq z(x, \omega)\leq z_{\max}$ with $z_{\min}$ and $z_{\max}$ representing two constants. Therefore (2.3) holds under this assumption. Furthermore, we assume $z(x, \omega)$ is a $\mathcal{F}_0$-measurable, mean-zero, Whittle-Matérn Gaussian random field, which is a stationary random field with the covariance function

where $\Gamma(\cdot)$ is the Gamma function, and ${g}(\cdot)$ stands for the inverse Fourier transform of $\hat g(\xi): = \frac{2^{q-\frac{1}{2}}\Gamma(q+\frac{1}{2})}{(1+\xi^2)^{q+\frac{1}{2}}}$. It is known that the parameter $q$ shown in (2.4) controls the regularity of the random field $z(x, \omega)$ ^[53]. Therefore, by taking different $q$ value, it's easy to numerically test and compare different cases. Notice that the covariance function and mean function uniquely determine a Gaussian random field ^[53].

Remark 2.1. The problem $(2.2)$ is a combination of an uncertainty quantification (UQ) problem with random diffusion coefficient field ^[47,68] and a stochastic model with multiplicative noise, which is a new problem in the sense that it involves two types of noise. We believe that such a problem can be applied to numerically investigate the effect of different types of perturbations on many physical models ^[47,48]. The log-Gaussian random field has received a lot of attention in the study of UQ problems ^[3,53], and appeared in some applications, e.g., geostatistical modelling ^[38,62]. The uniform boundedness condition imposed on the random variable $z(x, \omega)$ guarantees that the constants produced in the subsequent analysis is independent of $\omega$. We note that the so called trajectory rejection method proposed by Milstein et al. ^[55] uses similar idea of "truncation".

The theoretical result established in this paper depends on the following assumption on the nonlinear drift term $f(\cdot)$:

These assumptions are often used to prove the existence and uniqueness of the solution for SPDEs, see, e.g., ^[44,53].

We are interested in the mild solution of problem (2.2) in the Itô sense ^[17], defined by

where $S(t): = e^{-tL}$ is a semigroup generated by the operator $L$ ^[25]. The well-posedness of the problem (2.2) thus consists in verifying that the integrals in (2.7) are well defined and a function $u$ satisfying the integral equation (2.7) uniquely exists. We first notice that the realization of the random field $a(x, \omega)$ given in (2.3) is $2$ times mean-square differentiable due to $q > 2$ ^[53]. Thus, almost surely ($\mathbb{P}$-a.s.), $a(x, \omega)\in C^2(D)$. One verifies readily that $\mathcal{D}(L) = H^2(D)\cap H_0^1(D)$ almost surely ^[4], where $\mathcal{D}(L)$ is the domain of the operator $L$.

We also need some assumptions on the nonlinear term $G$, which are collected below:

- $L^{s}G(\cdot)$, $0\leq s\leq\frac{1}{2}$, is a mapping from $L^2(D)$ to ${\mathcal L}_Q$ such that:

- $\{G(v(\tau)):\tau\in[0, T]\}$ is a predictable ${\mathcal L}_Q$-valued process, such that

Remark 2.2. The assumptions $(2.8)$ and $(2.9)$ impose some restrictive conditions on the nonlinear term $G(\cdot)$, which include a combination of the nonlinear term $G(\cdot)$, the elliptic operator $L$, and the covariance operator $Q$. Notice that the similar or more general assumptions have been considered in ^[2,31,69].

We define the space $\mathbb{L}_2^t$ for $t\in [0, T]$, which is the Banach space of $L^2(D)$-valued predictable processes $\{v(\tau):\tau\in[0, t]\}$, equipped with the norm

Now we are in a position to state and prove the well-posedness of the mild solution to (2.2).

Theorem 2.1. Suppose that the initial value $u_0\in L^2(\Omega, L^2(D))$ is an $\mathcal{F}_0$-measurable random variable. Then, there exists a unique mild solution $u\in \mathbb{L}_2^T$ to $(2.2)$. Furthermore, the following stability inequality holds

Proof. We define the integral operator $\mathcal{M}$ by: for all $v\in \mathbb{L}_2^t$, $0\leq t\leq T$,

Obviously if there is a fixed point $u\in \mathbb{L}_2^T$ for the operator $\mathcal{M}$, then this fixed point is a mild solution defined by (2.7). The proof basically consists of two steps: 1) prove that the integral operator $\mathcal{M}$ is well-defined under the assumptions given above; 2) use the Fixed Point Theorem [53, Theorem 1.10] to establish the existence of a unique mild solution. This can be done by following the same lines as in [53, Theorem 10.26], using the imposed assumptions and a number of known results including the Karhunen-Loève (K-L) expansion of $Q$-Wiener process $W(s)$, It$\hat{\mathrm{o}}$ isometry, and the semigroup inequality [53, Exercise 10.7]: for $w\in L^2(\Omega, L^2(D)), \tau\in(0, T)$,

We emphasize here that $S(\tau)$ involves the random diffusion coefficient, thus the inequality (2.13) must be understood in the sense of almost surely. This, compared to the case of deterministic diffusion coefficient (see, e.g., [53, Theorem 10.26] for details), causes no essential difficulty in establishing the desired results.

3.   Random field sampling and fully discrete scheme

Our first goal in this section is to employ a method called approximation circulant embedding with padding ^[67] to uniformly sample the random diffusion coefficient $a(x, \omega)$. It is notable that some other sampling methods are available, such as turning bands method ^[22,54] and quadrature sampling method ^[60,61]. However the turning bands method is only applicable to isotropic Gaussian random fields, and the quadrature sampling method needs to know the spectral density function of the covariance function of random fields. One of the advantages of the sampling method we employ here is its applicability to stationary Gaussian random fields including isotropic random fields, and does not require prior knowledge of the spectral density function of the covariance function.

It is obvious from (2.3) that if we want to sample $a(x, \omega)$, we only need to sample $z(x, \omega)$. Notice that when $d = 1$, the resulting covariance matrix is Toeplitz, when $d = 2$ it is block Toeplitz, and when $d = 3$, it is nested block Toeplitz. The crucial ingredient of the circulant embedding sampling is that the target covariance matrix can be embedded into a large circulant matrix or a (nest) block circulant matrix, which can be decomposed by discrete Fourier transform. Then a new random field based on the combination of decomposition factors is constructed, which will be used to obtain the approximations of $z(x, \omega)$ for $x\in\bar{D}$. To be specificity and convenience, we briefly describe this approach by taking one-dimensional sampling as an example in this section.

Another purpose in the section is to present semi-implicit Euler-Maruyama scheme and finite element method to discrete problem (2.2) in time and space, respectively. We start by random field sampling.

3.1.   One-dimensional approximation circulant embedding with padding

Consider uniform sampling of random field $z(x, \omega)$ for $x\in[0, 1]$. Let $P$ be a positive integer. Set

Let $C: = (c_{ij})$ denote the $P\times P$ covariance matrix with respect to $z(x_p, \omega)$ for $p = 1, ..., P$, where $c_{ij}: = {\rm{cov}}(z(x_{i}, \omega), z(x_{j}, \omega)) = c_{q}(|x_{i}-x_{j}|)$ for $i, j = 1, ..., P$. Let $c_{i-j}: = c_{ij}$, then

One verifies readily that ${C}$ is a symmetric Toeplitz matrix, and it can be well defined by its first column $\boldsymbol{c_1} = (c_0, ..., c_{P-1})^T\in\mathbb{R}^P$. If we define $\boldsymbol{\bar{c}_1}: = \tbinom{\boldsymbol{c_1}}{\boldsymbol{0}}\in\mathbb{R}^{P+M}$ with $\boldsymbol{0}\in\mathbb{R}^{M}$ be a zero padding vector, a new symmetric Toeplitz matrix denoted by $\bar{C}\in \mathbb{R}^{(P+M)\times(P+M)}$ can be generated from $\boldsymbol{\bar{c}_1}$. Next, we carry out the minimal circulant extension [53, Definition 6.48] to $\bar{C}$ such that it can be embeded into a bigger circulant matrix denoted by $\tilde{\bar{{C}}}\in\mathbb{R}^{2\tilde{P}\times2\tilde{P}}$ for $\tilde{P}: = P+M-1$. Let $\tilde{\boldsymbol{\bar{c}}}_1$ be the first column of $\tilde{\bar{{C}}}$, $W^{*}$ represent the conjugate transpose of discrete Fourier matrix $W\in\mathbb{C}^{2\tilde{P}\times2\tilde{P}}$, and ${{d}}_j$ be the $j$-th entry of $\sqrt{2\tilde{P}}{{W}}^{*}\tilde{\boldsymbol{\bar{c}}}_1$. Then by Fourier representation, the circulant matrix $\tilde{\bar{{C}}}$ can be decomposed as follows:

where ${{\Lambda}}_{\pm}$ represents the diagonal matrix whose $j$-th diagonal element is $\pm\lambda_j: = \max\{0, \pm {{d}}_j\}$, i.e.,

Let $\boldsymbol{z}: = \big(z(x_1, \omega), \dots, z(x_{P}, \omega)\big)^T$. Our main goal is to take the sample approximations to the random vector $\boldsymbol{z}$. To this end, we construct a new random field vector $\boldsymbol{{Z}}$, defined by

where ${\rm{CN}}(\cdot, \cdot)$ denotes the complex Gaussian distribution [53, Definition 6.15]. It's readily to deduce that $\boldsymbol{{Z}}\sim{\rm{CN}}(\boldsymbol{0}, 2(\tilde{\bar{{C}}}+\tilde{\bar{{C}}}_{-}))$ with $\tilde{\bar{{C}}}_{-}: = {{W}}{{\Lambda}}_{-}{{W}}^{*}$, which means both real and imaginary parts of $\boldsymbol{{Z}}$ obey real Gaussian distribution N$(\boldsymbol{0}, (\tilde{\bar{{C}}}+\tilde{\bar{{C}}}_{-}))$. Notice that $\Vert \tilde{\bar{{C}}}_{-}\Vert_2\leq\rho({{\Lambda}}_{-})$ with $\rho({{\Lambda}}_{-})$ representing the spectral radius of ${{\Lambda}}_{-}$, and it is known that $\rho({{\Lambda}}_{-})$ can be small enough by increasing the dimension $M$ of zero padding vector ^[67]. Therefore $\tilde{\bar{{C}}}$ can be approximately treated as a non-negative definite matrix when the dimension $M$ is large enough, which is crucial for obtaining a good approximation of the random vector $\boldsymbol{z}$. Then the sample approximations of the random vector $\boldsymbol{z}$ can be provided by truncating the real or imaginary part of $\boldsymbol{Z}$.

The sampling procedure is summarized as follows:

i) Embed $C$ shown in (3.1) into the padded circulant matrix $\tilde{\bar{{C}}}\in\mathbb{R}^{2(P+M-1)\times2(P+M-1)}$ with dimension $M$ large enough;

ii) Compute $\Lambda_{+}$ by (3.2);

iii) Construct a new random field vector $\boldsymbol{Z}$ by (3.3) and take its real or imaginary part, denoted by $\boldsymbol{Z_1}\in\mathbb{R}^{2(P+M-1)}$;

iv) Truncate the first $P$ terms of $\boldsymbol{Z_1}$ and use it as an approximation to the random vector $\boldsymbol{z}$.

It is worthwhile to point out that the sampling method described above is convenient in the sense that it can simultaneously produce two sets of independent and identically distributed (i.i.d) samples in one sampling.

The sampling description for $d = 2$ or $d = 3$ is omitted here, and we refer to ^[67] for details. Notably, for each of the sampling data of the random diffusion coefficient, the problem (2.2) becomes a SEE with randomness only on the $G$-term.

3.2.   Spatio-temporal discretization

In this subsection we propose and analyze a discretization method for the problem (2.2). The proposed method is based on a finite element discretization in space and semi-implicit Euler-Maruyama approach in time.

We first briefly describe the $\mathbb{P}_1$ finite element method for spatial discretization. Let $\mathcal{T}_h$ be a shape-regular uniform partition of the domain $D$ such that $D = \cup_{K\in\mathcal{T}_h}K$, $K$ are intervals in the case $d = 1$, triangles in the case $d = 2$ and tetrahedrons for $d = 3$. Let $h_K$ be the diameter of the element $K$ and $h$ be the maximum of $\{h_K\}_{K\in\mathcal{T}_h}$. Assume that the partition $\mathcal{T}_h$ also equips the following characteristics:

● For different dimension $d$, the intersection of two different elements of $\mathcal{T}_h$, if not empty, is a vertex, a whole edge or a whole face of both of them;

● For $d = 2$ or $d = 3$, the ratio of the diameter $h_K$ of any element $K$ to the diameter of its inscribed circle ($d = 2$) or ball ($d = 3$) is smaller than a constant independent of $h$.

Define the finite element space $V_h$ by

where $\mathbb{P}_1(K)$ denotes the space of the polynomials of degree $\le$ 1 defined in $K$. Let $\mathcal{P}_{h}$ be the orthogonal projection from $L^2(D)$ to $V_h$, and $\mathcal{P}^{w}_{J}$ be the projection from $H^{\gamma}_0(D)$ to the finite-dimensional space span$\{\phi_1, \dots, \phi{_{J}}\}$. The spatial semi-discrete scheme of the problem (2.2) reads: find finite element approximation $u_h(t)\in V_h$ such that

where $L_h$: $V_h\to V_h$ is the finite-dimensional operator defined by

with $(\cdot, \cdot)$ be the $L^2(D)$-inner product.

We now describe the temporal discretization. Let $N$ be a positive integer, $\Delta t: = T/N$ be the uniform time step. Then the spatio-temporal full discretization of the problem (2.2), called hereafter the finite element method/semi-implicit Euler Maruyama scheme, reads:

where $\mathcal{P}^{w}_{J}\Delta W^n: = \sum_{j = 1}^{J} \sqrt{q_j}(\beta_j(t_{n+1})-\beta_j(t_{n})) \phi_j$ with $\beta_j(t)$ be the i.i.d $\mathcal{F}_t$-Brownian motions.

In actual calculation, we will use the average of the sampled values at the nodes of the element $K$ to approximate $a(x, \omega)$ for $x\in K$. Therefore, the overall cost of the scheme (3.5) is roughly equal to solving a linear equation with random variable coefficients of the form $(I+{\Delta t} L_h)\boldsymbol{x} = \boldsymbol{b}$ at each time step.

4.   Error estimate

This section is devoted to analyzing the strong convergence error of the spatio-temporal full discretization (3.5) to the mild solution (2.7). Here, strong convergence is understood in the sense of convergence with respect to the norm $\|\cdot\|_{L^2(\Omega, L^2(D))}$. We first note that the full-discrete scheme (3.5) can be rewritten under form:

It is readily seen that $L_h$ is reversible in $V_h$, i.e., $L_h^{-1}v_h$ is well defined for all $v_h\in V_h$. We now extend the definition of $L_h^{-1}$ to all $v\in L^2(D)$ by $L_h^{-1}v = L_h^{-1}\mathcal{P}_{h}v$. By the assumption on $a(x, \omega)$, we know that for almost every $\omega\in\Omega$, $L_h^{-1}$ is a non-negative definite operator from $L^2(D)$ to $V_h$. In fact, for all $v\in L^2(D)$, there exists $w_h\in V_h$ such that $L_h w_h = \mathcal{P}_{h}v$, and thus

Let $S_{h, \Delta t}^{n}: = (I+\Delta tL_h)^{-n}$. The fully discrete approximation can be expressed under the form:

Subtracting (4.2) from the mild solution (2.7) gives

with $\theta_i, i = 1, 2, 3$, representing

Our goal in the following is to estimate $\theta_1$, $\theta_2$, $\theta_3$ separately in the sense of strong convergence. To this end, we first give some preliminaries that will be used in subsequent analysis.

$\bullet$ If the initial value $u_0\in L^2(\Omega, \mathcal{D}(L))$ is an $\mathcal{F}_0$-measurable random variable, then there exists a constant $c$ depended on $u_0$ such that the mild solution $u$ defined in (2.7) satisfies the following temporal Hölder regularity:

The proof of (4.7) can be done by following the same lines as in [53, Lemma 10.27], which is omitted here. Basically, it makes use of the properties of the operator $L$ and its associated semigroup $S(t)$, satisfied in the sense of almost surely.

$\bullet$ The operator $L$ and the induced semigroup $S(t)$ satisfy the following estimates, which is a straightforward extension of the classical results (see, e.g., ^[44,64]) to the sense of almost surely:

- For each $\alpha\ge 0$, there exists a constant $c$ such that

- For $\alpha\in[0, 1]$, there exists a constant $c$ such that

$\bullet$ The nonlinear term $G$ satisfies

In fact, it follows from the assumption (2.9): for all $v_1, v_2\in{L^2(D)}$,

4.1.   Strong error estimate

We first focus on strong error estimate for the term $\theta_1$. It is worth pointing out that, although our analysis is inspired by the work ^[44] based on the rational function approach, our proof makes full use of the standard framework of the finite element approximation to the linear parabolic equation as well as the fact that operator $L_h^{-1}$ is non-negative definite from $L^2(D)$ to $V_h$. Let

Then it follows from the definition (4.4) that $\theta_1 = T_nu_0$.

Lemma 4.1 (Error estimate of $\theta_1$). Suppose $u_0\in L^2(\Omega, \mathcal{D}(L))$ is an $\mathcal{F}_0$-measurable random variable. Then there exists a constant $c$ independent of $h$ and $\Delta t$ (but depends on $u_0$), such that

where $\theta_1$ is given in (4.4).

Proof. Obviously, $\theta_1$ characterizes the error between the exact solution $e^{-t_nL}u_0$ and the full discrete solution $(I+\Delta tL_h)^{-n}\mathcal{P}_{h}u_0$, which can be splited into two parts:

where

is the spatial discretization error, while

is the temporal discretization error. Clearly, we have $e^n_1 = y(t_n) - y_h(t_n)$, where $y(t)\in H_0^1(D)$ and $y_h(t) \in V_h$ are the solutions of the parabolic equation

and its finite element semi-discrete equation

respectively. Let

It can be verified that

Using the non-negative definite of the operator $L_h^{-1}$ as well as the standard error analysis of the finite element approximation to the parabolic equation [53, Lemma 3.51] gives: for almost every $\omega\in\Omega$,

According to $y(\tau) = e^{-L\tau}u_0$, and the standard analytical framework for the finite element approximation to the linear parabolic equation, the terms in the right-hand side can be bounded by:

Thus

We now turn to estimate the temporal discretization error $e^n_2$. A direct calculation gives

Noticing that the operator $L_h$ is symmetric, and the ${\mathcal L}(L^2(D))$-norm of the operator $\big(e^{-n\Delta tL_h}-(I+\Delta tL_h)^{-n}\big)L_h^{-1}$ is equal to its spectral radius, i.e.,

where $\lambda^h_{j} > 0$, $j = 1, ..., K$, are the eigenvalues of $L_h$. Note that $|(e^{-nx}-(1+x)^{-n})/x|$ is bounded for $x > 0$, therefore taking $x = \lambda_{j}^h\Delta t$ gives

This proves

Combining (4.12), (4.13), and (4.14) gives

The above estimate holds for almost all $\omega\in \Omega$. Therefore

Remark 4.1. If the $\mathcal{F}_0$-measurable random variable $u_0\in L^2(\Omega, L^2(D))$. Then we have only ^[44,53]: for almost every $\omega\in\Omega$,

We next derive the error estimate for the term $\theta_2$, which is based on the standard error analysis for the deterministic semilinear evolution equation, the semigroup property, and the temporal regularity of the mild solution.

Lemma 4.2 (Error estimate of $\theta_2$). Suppose $u_0\in L^2(\Omega, \mathcal{D}(L))$ is an $\mathcal{F}_0$-measurable random variable. Then there exists a constant $c$ independent of $h$ and ${\Delta t}$ (but depends on $\|u_0\|_{L^2(\Omega, \mathcal{D}(L))}$), such that

where $\theta_2$ is given by $(4.5)$.

Proof. The term to be bounded can be expressed by

which can be decomposed into

with

For the part $\theta^1_2$, it follows from the norm definition (2.1):

According to (4.8) and (4.9), the operator norm $\| S(t_n-\tau)-S(t_n-t_k)\|_{_{\mathcal{L}({L^2(D)})}}$ is bounded $\mathbb{P}$-a.s. by:

We further use (2.5) and (2.11) to derive

The estimate of $\theta^2_2$ follows from (4.11), (4.15), (2.5), and (2.11):

For the part $\theta^3_2$, by $\| S_{h, \Delta t}^{n-k}\|_{_{\mathcal{L}({L^2(D)})}}\leq1$ ($\mathbb{P}\mbox{-a.s.}$), $\| \mathcal{P}_{h}\|_{_{\mathcal{L}({L^2(D)})}}\leq1$, (2.6), and (4.7), we have

The part $\theta^4_2$ can be estimated similarly:

Finally, we conclude by combining all above estimates with the triangle inequality.

In order to estimate the error contribution term $\theta_3$, we need to derive an estimate related to the nonlinear term $G$.

Lemma 4.3. Suppose $u_0\in L^2(\Omega, \mathcal{D}(L))$ is a $\mathcal{F}_0$-measurable random variable, and the eigenvalues of $Q$ satisfy $q_j = \mathcal{O}(j^{-(2\gamma+1+{\epsilon})})$ for some $\gamma\ge0$ and ${\epsilon} > 0$. Then it holds: for $0\leq \tau_1\leq \tau_2\leq T$,

Proof. Using the triangle inequality:

we are led to estimate the two terms in the right-hand side. First, employing (2.9) and (4.7) gives:

Then under the assumptions (2.8) and (2.11), we have

This proves (4.16).

Lemma 4.4 (Error estimate of $\theta_3$). Under the assumptions of Lemma 4.3, further assume ${\Delta t} = O(h^2) = O(J^{-\gamma})$. Then there exists a constant $c$ independent of ${\Delta t}$ and $h$, such that

where $\theta_3$ is given by $(4.6)$.

Proof. Split $\theta_3$ as $\theta_3 = \sum_{i = 1}^{4}\theta_3^i$, where

with

For the part $\theta_3^{1}$, it follows from the Itô isometry:

We are led to estimate the two terms on the right-hand side of the inequality. First using $\big\|S(t_n-\tau)-S(t_n-t_{n-1})\big\|^2_{_{{\mathcal L}(L^2(D))}}\leq c$ ($\mathbb{P}$-a.s.), (2.8), and (2.11) yields

Then employing (4.8) and (4.9) gives

Making use of (2.8), (2.11) gives

Therefore

The estimate for $\theta_3^{2}$ follows from Itô isometry, (4.11), and (4.15):

We further use (2.8), (2.11), and ${\Delta t} = O(h^2)$ to get

For the part $\theta_3^{3}$, employing $O(h^2) = O(J^{-\gamma})$, $\| S^{n-k}_{h, \Delta t}\|_{_{\mathcal{L}({L^2(D)})}}\leq1$ ($\mathbb{P}$-a.s.), and (4.16) gives

The last part $\theta_3^{4}$ can be estimated by using Itô isometry, $\| S^{n-k}_{h, \Delta t}\|_{_{\mathcal{L}({L^2(D)})}}\leq1$ ($\mathbb{P}$-a.s.), and (4.10):

Finally we combine all above estimates and keep only the leading order to conclude.

Thanks to the results established in the previous lemmas, we are now in a position to derive the full discretization error bound, which is stated in the following theorem.

Theorem 4.1. Let $u$ be the mild solution defined in $(2.7)$, and $u_h^n$ be the numerical solution of $(3.5)$. Then under the assumptions stated in Lemmas 4.1–4.4, there exists a constant $c$ independent of ${\Delta t}$ and $h$, such that

Proof. It follows from (4.3), Lemmas 4.1-4.4, and the triangle inequality:

Then the discrete Gronwall inequality yields

This ends the proof.

Remark 4.2. Notice that the term $\Delta t^{\frac{1}{2}}$ dominates the term $\Delta t \ln(\Delta t^{-1})$, the estimate given in Theorem 4.1 can be simplified by

Also notice that $\Delta t^{-\varepsilon_0}$ dominates $\ln(\Delta t^{-1})$ for arbitrarily small $\varepsilon_0 > 0$, we have

or, since $\Delta t = O(h^2)$,

5.   Numerical results

Several numerical experiments are presented in this section to validate our theoretical estimates and show the effect of stochastic factors on numerical solutions. We start by testing the convergence orders of time and space.

Example 5.1 (Accuracy Test). We take the one-dimensional stochastic Langmuir reaction equation with random diffusion coefficient field and multiplicative force noise as a numerical example to test the temporal and spatial convergence orders. The underlying equation is expressed as:

where $z(x, \omega)$ is a bounded Gaussian random field with mean-zero and covariance function $c_q(x)$, and $W(x, t)$ is a $H^{\gamma}_0$-valued Wiener process defined by

where $q_j = \mathcal{O}(j^{-(2\gamma+1+\epsilon)})$ with arbitrary small positive $\epsilon$.

We first test the effectiveness of the sampling method used in this work. It is known from ^[67] that if the spectral radius $\rho({{\Lambda}}_{-})$ of the diagonal matrix ${\Lambda}_{-}$ defined in (3.2) is bigger than zero, then it will rapidly decrease to zero as the dimension $M$ of the zero padding vector increases. We run the sampling method presented in Subsection 3.1 to produce an approximation of the random vector $\boldsymbol{z} = \big(z(x_1, \omega), \dots, z(x_{P}, \omega)\big)^T$ by taking $P = 100$ and $q = 3$. Table 1 shows the spectral radius of ${\Lambda}_{-}$ as a function of $M$. As expected, we observe that $\rho({{\Lambda}}_{-})$ rapidly decreases to a level close to the machine precision as $M$ increases. Furthermore, taking $M = 4000$, we calculate the mean values of the sampled values of $\boldsymbol{z}$ for different sample numbers. The obtained result is given in Figure 1, from which we observe that as the sampling number increases, the mean value of the sampling approximation of the random field $\boldsymbol{z}$ will gradually approach the theoretical mean-zero function. For each $x\in\bar{D}$, we denote by $\mu_m(x)$ the approximation of $\mathbb{E}[z(x, \omega)]$ under $m$ samples. We calculate $\mu_m(x)$ for $x = 0.1, 0.5, 0.9$ by taking $q = 3, M = 4000$. The obtained results for $m = 10, 20, 50, 100$ are shown in Figure 1, from which we observe that $|\mu_m(x)|$ converges to the theoretical mean-zero as the sampling number increases, and convergence rate is roughly $O(m^{-\frac{1}{2}})$.

The strong convergence rate in space and time is measured in terms of mean-square approximation errors at the endpoint $T = 0.1$, caused by the spatial and temporal discretizations. The expected value of error is approximated by computing the mean of 200 samples. Note that the exact solution of the problem (5.1) is unknown, and we will use the reference solution computed in the fine space-time mesh size as an approximation to the exact solution. If we denote by $u^{\rm{ref}}_j$ the reference solution of the $j$-th sample of the exact solution $u(T)$, and denote by $u_{j, h}^N$ the value of the $j$-th sample of the fully discrete numerical solution $u_h^N$. Then the mean-square error $\|u(T)-u^N_h\|_{_{L^2(\Omega, L^2(D))}}$ is approximately calculated by

We first test the time accuracy with different nonlinear terms $G$. Take the numerical solution computed by spatial mesh $h = 1/128$ and time step size $\Delta t = 10^{-6}$ as the reference solution for every sample. The approximation error $u_{\rm{error}}$ under different time step size is calculated by taking $u_0(x) = \sin(2\pi x)$, $\varepsilon = 10^{-3}$, $\gamma = 1$ and $q = 2$. Tables 2 and 3 respectively show the results for the cases where $G(u) = (1-u^2)/2$ and $G(u) = u/2$, from which we observe that this is as predicted by the theory.

Next, we test the spatial accuracy. Now take the numerical solution computed by $h = 1/512$ and $\Delta t = 10^{-5}$ as the reference solution for every sample. We compute the approximation error $u_{\rm{error}}$ under different spatial mesh size by taking $u_0(x) = \sin(2\pi x)$, $\varepsilon = 10^{-3}$, $\gamma = 1$ and $q = 2$ again. Tables 4 and 5 separately shows the relevant error and spatial convergence order for $G(u) = (1-u^2)/2$ and $G(u) = u/2$, which is also consistent with the theoretical result.

Example 5.2 (One-dimensional phenomenon comparison). In this example, the time evolution of the numerical solution of the stochastic Allen-Cahn equation is compared to that of the deterministic Allen-Cahn equation to show the effect of random perturbations, where the deterministic version is expressed by:

We first show the effect of random field $a(x, w) = \varepsilon e^{z(x, w)}$ on the numerical solution in the absence of the nonlinear term $G$ (i.e., $G(u) = 0$), where $z(x, \omega)$ is a bounded mean-zero Gaussian random field with covariance function $c_q(x)$. Given a sample point, by taking $u_0 = \sin(4 \pi x)$, $T = 0.1$, $\Delta t = 10^{-5}$, $h = 1/128$ and $\varepsilon = 10^{-2}$, we plot in Figure 2 the contour maps of the numerical solution under different cases, where figure (a) represents the deterministic case and figures (b) and (c) denote the random case with $q = 0.1$ and $q = 2$, respectively. Compared to the deterministic model, it is seen from figures (b) and (c) that the random diffusion coefficient makes the diffusion process uncertain. Moreover, it's known that the larger the parameter $q$, the more regular the random field $z(x, \omega)$ ^[53], which results in the diffusion process shown in figure (c) is more uniform than that in figure (b).

Then we give a demonstration of the case with both random diffusion coefficients as well as multiplicative force noise. Given a sample point, by taking $u_0 = \sin(4 \pi x)$, $T = 4$, $\Delta t = 10^{-4}$, $h = 1/128$, $q = 2$, $\varepsilon = 10^{-5}$ and $G(u) = \frac{1}{2}(1-u^2)$, we plot the time evolution of the numerical solution in Figure 3, where figure (a) denotes the deterministic model and figures (b) and (c) represent the case with $\gamma = 0.5$ and $\gamma = 1$, respectively. Compared to figure (a), it can be seen from figures (b) and (c) that there are small-scale structures resulted from noise, which are not present in the deterministic model. Notably, the static kink corresponding to the deterministic model varies greatly after the incorporation of noise and random diffusion coefficient fields. The kinks can interact, even annihilate each other, and some new kinks may arise. One more thing to point out that the larger the regularity parameter $\gamma$, the smoother the noise and the smaller the kink variation, which seems to be observed between figures (b) and (c).

Example 5.3 (Two-dimensional phenomenon comparison). In this example, numerical simulations of surfaces evolving according to the mean curvature of two-dimensional deterministic and stochastic Allen-Cahn (AC) equation are presented to show the perturbing effects of random factors on the numerical solution. The underlying equations of deterministic and random versions are shown as below:

where $z(x, \omega)$ stands for a bounded two-dimensional Gaussian random field with mean-zero and covariance function $c_q(x)$. Let $x = (x_1, x_2)^T$, $x_1, x_2\in(0, 1)$. Denote $W(x, t)$ by

where $q_{j_1, j_2} = \exp(-\frac{j_1^2+j_2^2}{200})$ and $\beta_{j_1, j_2}(t)$ are the i.i.d Brownian motions. The equations shown in (5.2) are supplemented with the initial condition $u(x, 0) = u_0(x)$, and the homogeneous Neumann boundary condition.

For the deterministic AC equation, it is known that as $\varepsilon \rightarrow 0$, the zero level set of $u$, which is denoted by $\Gamma_{t}^{\varepsilon}: = \{x\in D:u(x, t) = 0\}$, approaches a surface $\Gamma_{t}$ whose evolution follows the geometric law:

where $V$ is the normal velocity of the surface $\Gamma_{t}$ at each point, $\kappa$ is its mean curvature, and $R$ is the principal radii of curvature at the point of the surface ^[46,49]. If we denote the radius at time $t$ by $R(t)$ and set the initial radii to be $R_0$, then $R(t) = \sqrt{R_{0}^{2}-2 t}$ due to $V(t) = \frac{d R(t)}{d t} = -\frac{1}{R(t)}$.

The first row of the Figure 4 shows the evolution of a star-shaped interface in a curvature-driven flow of the deterministic AC model on the computational domain $D = (0, 1)\times(0, 1)$ with a 512$\times$512 mesh and $\varepsilon = 7.5\times10^{-4}$, ${\Delta t} = 5\times 10^{-5}$. The initial configuration is defined by:

It is observed that the tips of the star move inward, while the gaps between the tips move outward, and the whole shape shows a trend of shrinking towards the center.

Given a sample point, the second row of Figure 4 shows the effect of the random diffusion coefficient field on the numerical solution for the case of $q = 0.5$ and $G(u) = 0$. We see that the evolution of the initial star interface is no longer shrink regularly, and its shape changes randomly.

Finally, for a fixed sample point, we show on the third row of Figure 4 the evolution of the star interface for the case with $q = 0.5$ and $G(u) = 10(1-|u|)$. It is seen that noise plays a significant role, it changes the properties of the solutions. Similar to the case in one-dimensional, there are small-scale structures generated by noise. The kinks can interact or even annihilate each other and some new kinks appear.

6.   Conclusions

In this paper we study a stochastic evolution equation driven by a bounded $\log$-Whittle-Mat$\acute{{\mathrm{e}}}$rn random diffusion coefficient field and $Q$-Wiener multiplicative force noise. The well-posedness of the underlying equations is established, and an efficient numerical method is proposed and analyzed. The proposed method makes use of a circulant embedding approach with padding to sample the random coefficient field, and semi-implicit Euler-Maruyama scheme and finite element method for the spatio-temporal discretization. An estimate for the strong convergence rate is derived. Several numerical examples are provided to validate the proposed method.

Acknowledgments

This research is partially supported by NSFC grants 11971408 and 12271457, NSFC/ANR joint program ANR-16-CE40-0026-01, and the French State in the frame of the "Investments for the future" programme Idex Bordeaux ANR-10-IDEX-03-02.