<figure><img src="image_1.png" /></figure>

HAL
open science

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

Xiao Qi, Mejdi Azaiez, Can Huang, Chuanju Xu

► To cite this version:

Xiao Qi, Mejdi Azaiez, Can Huang, Chuanju Xu. An efficient numerical approach for stochastic evolution PDEs driven by random diffusion coefficients and multiplicative noise. AIMS Mathematics, 2022, 7 (12), pp.20684-20710. 10.3934/math.20221134. hal-03962943v2

HAL Id: hal-03962943

https://hal.science/hal-03962943v2

Submitted on 3 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Research article

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

Xiao Qi<fnref n="1" />, Mejdi Azaiez<fnref n="1" />,<fnref n="2" />, Can Huang<fnref n="1" /> and Chuanju Xu<fnref n="1" />,.*

* Correspondence: Email: cjxu@xmu.edu.cn.

**Abstract:** In this paper, we investigate the stochastic evolution equations (SEEs) driven by a bounded log-Whittle-Matérn (W-M) random diffusion coefficient field and Q-Wiener multiplicative force noise. First, the well-posedness of the underlying equations is established by proving the existence, uniqueness, and stability of the mild solution. A sampling approach called approximation circulant embedding with padding is proposed to sample the random coefficient field. Then a spatio-temporal discretization method based on semi-implicit Euler-Maruyama scheme and finite element method is constructed and analyzed. An estimate for the strong convergence rate is derived. Numerical experiments are finally reported to confirm the theoretical result.

**Keywords:** SEEs; random coefficient; Q-Wiener multiplicative noise; strong convergence

**2010 Mathematics Subject Classification:** 60H15, 60H35, 65C50

**1. Introduction**

Stochastic partial differential equations (SPDEs) appear 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

<footnote n="1">¹ School of Mathematical Sciences and Fujian Provincial Key Laboratory of Mathematical Modeling and High Performance Scientific Computing, Xiamen University, 361005 Xiamen, China</footnote>

<footnote n="2">² Bordeaux INP, Laboratoire I2M UMR 5295, 33607 Pessac, France</footnote>
on the other hand, see, e.g., [7, 11, 13–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–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–15, 18, 24, 26, 29, 30,
34–36, 39, 43, 44, 51, 52, 59, 63–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á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:

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

* 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.

• 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 $O(h^{2-\epsilon_0}+\Delta t^{\frac{1}{2}})$ of the full discretization under certain assumptions, where $\epsilon_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 \ge 0$. $\mathcal{L}(L^2(D))$ represents the space of bounded linear operators $A: L^2(D) \to L^2(D)$ equipped with operator norm $\|A\|_{L(L^2(D))} = \sup_{u \neq 0} \frac{\|Au\|_{L^2(D)}}{\|u\|_{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 \ge 0}$. Additionally, we denote $v(x, \omega) \in L^2(\Omega, L^2(D))$ if

$$ \|v\|_{L^2(\Omega, L^2(D))} < +\infty, $$

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

$$ \|v\|_{L^2(\Omega, L^2(D))} := [\mathbb{E}\|v(\cdot, \omega)\|_{L^2(D)}^2]^{\frac{1}{2}} \quad (2.1) $$

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

$$ \left( \sum_{j=1}^{\infty} \| B Q^{\frac{1}{2}} \phi_j \|_{L^2(D)}^2 \right)^{\frac{1}{2}} < +\infty. $$

$\mathcal{L}_Q$ endowed with the norm $\|B\|_{\mathcal{L}_Q} := \left( \sum_{j=1}^{\infty} \| B Q^{\frac{1}{2}} \phi_j \|_{L^2(D)}^2 \right)^{\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

$$ \|B(\omega)\|_{L^2(\Omega, \mathcal{L}_Q)} := \mathbb{E}[\|B(\omega)\|_{\mathcal{L}_Q}^2]^{\frac{1}{2}}. $$
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:

$$ 
\begin{aligned}
du(x,t) &= (-Lu + f(u))dt + G(u)dW(x,t), \quad 0 < t < T, \quad x \in D, \\
u(x,t) &= 0, \quad 0 \le t \le T, \quad x \in \partial D, \\
u(x,0) &= u_0(x), \quad x \in \bar{D},
\end{aligned}
\quad (2.2) $$

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 $\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$,

$$ 0 < a_{\min} \le a(x, \omega) = \varepsilon e^{z(x, \omega)} \le a_{\max} < \infty. \quad (2.3) $$

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} \le z(x, \omega) \le 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 $F_0$-measurable, mean-zero, Whittle-Matérn Gaussian random field, which is a stationary random field with the covariance function

$$ c_q(\|x\|_2) := \frac{g(\|x\|_2)}{2^{q-1}\Gamma(q)}, \quad x \in \bar{D}, \quad q > 2, \quad (2.4) $$

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)$:

$$ \|f(v)\|_{L^2(D)} \le c(1 + \|v\|_{L^2(D)}), \quad \forall v \in L^2(D), \quad (2.5) $$

$$ \|f(v_1) - f(v_2)\|_{L^2(D)} \le c(\|v_1 - v_2\|_{L^2(D)}), \quad \forall v_1, v_2 \in L^2(D). \quad (2.6) $$
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

$$u(t) = S(t)u_0 + \int_0^t S(t-\tau)f(u(\tau))d\tau + \int_0^t S(t-\tau)G(u(\tau))dW(\tau), \quad (2.7)$$

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 (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 \le s \le \frac{1}{2}$, is a mapping from $L^2(D)$ to $\mathcal{L}_Q$ such that:

$$\|L^s G(v)\|_{\mathcal{L}_Q} \le c(1 + \|v\|_{L^2(D)}), \quad \forall v \in L^2(D), \quad (2.8)$$

$$\|L^s(G(v_1) - G(v_2))\|_{\mathcal{L}_Q} \le c\|v_1 - v_2\|_{L^2(D)}, \quad \forall v_1, v_2 \in L^2(D). \quad (2.9)$$

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

$$\int_0^T \mathbb{E}[\|G(v)\|_{\mathcal{L}_Q}^2] d\tau < +\infty, \quad \forall v \in L^2(D). \quad (2.10)$$

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

$$\|v\|_{\mathbb{L}_2^t} := \sup_{\tau \in [0,t]} \|v(\tau)\|_{L^2(\Omega, L^2(D))} < +\infty.$$

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

$$\sup_{t \in [0, T]} \|u(t)\|_{L^2(\Omega, L^2(D))} \le c_T (1 + \|u_0\|_{L^2(\Omega, L^2(D))}). \quad (2.11)$$

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

$$(\mathcal{M}v)(t) := S(t)u_0 + \int_0^t S(t-\tau)f(v(\tau))d\tau + \int_0^t S(t-\tau)G(v(\tau))dW(\tau). \quad (2.12)$$
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ô isometry, and the semigroup inequality [53, Exercise 10.7]: for $w \in L^2(\Omega, L^2(D))$, $\tau \in (0, T)$,

$$ \|S(\tau)w\|_{L^2(\Omega, L^2(D))} \le \|e^{-\tau L}\|_{L^2(L^2(D))} \|w\|_{L^2(\Omega, L^2(D))} \le \|w\|_{L^2(\Omega, L^2(D))}. \quad (2.13) $$

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. $\square$

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

$$ 0 = x_1 \le \dots \le x_P = 1, \quad x_{P+1} - x_P = \frac{1}{P-1}, \quad p = 1, \dots, P-1. $$
Let $C := (c_{ij})$ denote the $P \times P$ covariance matrix with respect to $z(x_p, \omega)$ for $p = 1, \dots, P$, where $c_{ij} := \operatorname{cov}(z(x_i, \omega), z(x_j, \omega)) = c_q(|x_i - x_j|)$ for $i, j = 1, \dots, P$. Let $c_{i-j} := c_{ij}$, then

$$ C = \begin{pmatrix} c_0 & c_{-1} & \cdots & c_{1-P} \\ c_1 & c_0 & \cdots & c_{2-P} \\ \vdots & \vdots & \ddots & \vdots \\ c_{P-1} & \cdots & c_1 & c_0 \end{pmatrix}. \qquad (3.1) $$

One verifies readily that $C$ is a symmetric Toeplitz matrix, and it can be well defined by its first column $c_1 = (c_0, \dots, c_{P-1})^T \in \mathbb{R}^P$. If we define $\tilde{c}_1 := \begin{pmatrix} c_1 \\ 0 \end{pmatrix} \in \mathbb{R}^{P+M}$ with $0 \in \mathbb{R}^M$ be a zero padding vector, a new symmetric Toeplitz matrix denoted by $\tilde{C} \in \mathbb{R}^{(P+M)\times(P+M)}$ can be generated from $\tilde{c}_1$. Next, we carry out the minimal circulant extension [53, Definition 6.48] to $\tilde{C}$ such that it can be embedded into a bigger circulant matrix denoted by $\tilde{C} \in \mathbb{R}^{2\tilde{P}\times2\tilde{P}}$ for $\tilde{P} := P + M - 1$. Let $\tilde{c}_1$ be the first column of $\tilde{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{c}_1$. Then by Fourier representation, the circulant matrix $\tilde{C}$ can be decomposed as follows:

$$ \tilde{C} = W(\Lambda_+ - \Lambda_-)W^*, $$

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

$$ \Lambda_{\pm} = \operatorname{diag}(\pm\lambda_1, \dots, \pm\lambda_{2\tilde{P}}). \qquad (3.2) $$

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

$$ Z := W\Lambda_+^\frac{1}{2}\xi, \quad \xi \sim \operatorname{CN}(0, 2I_{2\tilde{P}}), \qquad (3.3) $$

where $\operatorname{CN}(\cdot, \cdot)$ denotes the complex Gaussian distribution [53, Definition 6.15]. It's readily to deduce that $Z \sim \operatorname{CN}(0, 2(\tilde{C} + \tilde{C}_-))$ with $\tilde{C}_- := W\Lambda_-W^*$, which means both real and imaginary parts of $Z$ obey real Gaussian distribution $\operatorname{N}(0, (\tilde{C} + \tilde{C}_-))$. Notice that $\|\tilde{C}_-\|_2 \le \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{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 $z$. Then the sample approximations of the random vector $z$ can be provided by truncating the real or imaginary part of $Z$.

The sampling procedure is summarized as follows:

i) Embed $C$ shown in (3.1) into the padded circulant matrix $\tilde{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 $Z$ by (3.3) and take its real or imaginary part, denoted by $Z_1 \in \mathbb{R}^{2(P+M-1)}$;

iv) Truncate the first $P$ terms of $Z_1$ and use it as an approximation to the random vector $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

$$V_h := \{v \in C^0(\bar{D}) \text{ with } v = 0 \text{ on } \partial D \text{ and } v|_K \in \mathbb{P}_1(K) \text{ for all } K \in \mathcal{T}_h\},$$

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}_h^w$ be the projection from $H_0^\gamma(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

$$du_h(t) = (-L_h u_h(t) + \mathcal{P}_h f(u_h(t)))dt + \mathcal{P}_h(G(u_h(t))\mathcal{P}_h^w dW(t)), \quad \forall 0 < t \le T, \\ u_h(0) = \mathcal{P}_h u_0, \qquad (3.4)$$

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

$$(L_h w, v) := (a(x, \omega) \nabla w, \nabla v), \quad \forall w, v \in V_h$$

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:

$$(I + \Delta t L_h)u_h^{n+1} = u_h^n + \Delta t \mathcal{P}_h f(u_h^n) + \mathcal{P}_h(G(u_h^n)\mathcal{P}_h^w \Delta W^n), \quad n = 0, \dots, N-1, \\ u_h^0 = \mathcal{P}_h u_0, \qquad (3.5)$$

where $\mathcal{P}_h^w \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)x = 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:

$$u_h^{n+1} = (I + \Delta t L_h)^{-1} (u_h^n + \Delta t \mathcal{P}_h f(u_h^n) + \mathcal{P}_h G(u_h^n) \mathcal{P}_j^w \Delta W^n), \quad n = 0, \dots, N-1. \quad (4.1)$$

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

$$\begin{aligned} (L_h^{-1}v, v) &= (L_h^{-1}\mathcal{P}_h v, v) = (L_h^{-1}\mathcal{P}_h v, \mathcal{P}_h v) = (L_h^{-1}L_h w_h, L_h w_h) \\ &= (w_h, L_h w_h) = (a(x, \omega)\nabla w_h, \nabla w_h) \ge 0. \end{aligned}$$

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

$$u_h^n = S_{h, \Delta t}^n \mathcal{P}_h u_0 + \sum_{k=0}^{n-1} \Delta t S_{h, \Delta t}^{n-k} \mathcal{P}_h f(u_h^k) + \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} S_{h, \Delta t}^{n-k} \mathcal{P}_h G(u_h^k) \mathcal{P}_j^w dW(\tau). \quad (4.2)$$

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

$$u(t_n) - u_h^n = \theta_1 + \theta_2 + \theta_3 \quad (4.3)$$

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

$$\theta_1 := S(t_n)u_0 - S_{h, \Delta t}^n \mathcal{P}_h u_0, \quad (4.4)$$

$$\theta_2 := \sum_{k=0}^{n-1} \left( \int_{t_k}^{t_{k+1}} S(t_n - \tau) f(u(\tau)) d\tau - \Delta t S_{h, \Delta t}^{n-k} \mathcal{P}_h f(u_h^k) \right), \quad (4.5)$$

$$\theta_3 := \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \left( S(t_n - \tau)G(u(\tau)) - S_{h, \Delta t}^{n-k} \mathcal{P}_h G(u_h^k) \mathcal{P}_j^w \right) dW(\tau). \quad (4.6)$$

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.

* 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:

$$\|u(\tau_2) - u(\tau_1)\|_{L^2(\Omega, L^2(D))} \le c(\tau_2 - \tau_1)^{\frac{1}{2}}, \quad \forall 0 \le \tau_1 \le \tau_2 \le T. \quad (4.7)$$

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.
• 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

$$ \|L^{\alpha} S(t)\|_{L^2(L^2(D))} \le ct^{-\alpha}, \quad \forall t > 0. \qquad (4.8) $$

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

$$ \|L^{-\alpha}(I - S(t))\|_{L^2(L^2(D))} \le ct^{\alpha}, \quad \forall t \ge 0. \qquad (4.9) $$

• The nonlinear term $G$ satisfies

$$ \|G(v_1)\mathcal{P}_J^w - G(v_2)\mathcal{P}_J^w\|_{L_Q} \le c\|v_1 - v_2\|_{L^2(D)}, \quad \forall v_1, v_2 \in L^2(D). \qquad (4.10) $$

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

$$ \begin{aligned} \|G(v_1)\mathcal{P}_J^w - G(v_2)\mathcal{P}_J^w\|_{L_Q} &= \left[ \sum_{j=1}^{\infty} \left\| (G(v_1) - G(v_2))\mathcal{P}_J^w Q^{1/2} \phi_j \right\|_{L^2(D)}^2 \right]^{\frac{1}{2}} \\ &= \left[ \sum_{j=1}^{J} \left\| (G(v_1) - G(v_2))Q^{1/2} \phi_j \right\|_{L^2(D)}^2 \right]^{\frac{1}{2}} \\ &\le \|G(v_1) - G(v_2)\|_{L_Q} \le c\|v_1 - v_2\|_{L^2(D)}. \end{aligned} $$

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

$$ T_n := (e^{-t_n L} - (I + \Delta t L_h)^{-n} \mathcal{P}_h). \qquad (4.11) $$

Then it follows from the definition (4.4) that $\theta_1 = T_n u_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

$$ \|\theta_1\|_{L^2(\Omega, L^2(D))} \le c(\Delta t + h^2), $$

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

*Proof.* Obviously, $\theta_1$ characterizes the error between the exact solution $e^{-t_n L} u_0$ and the full discrete solution $(I + \Delta t L_h)^{-n} \mathcal{P}_h u_0$, which can be split into two parts:

$$ \theta_1 = e_1^n + e_2^n, \qquad (4.12) $$

where

$$ e_1^n := e^{-t_n L} u_0 - e^{-t_n L_h} \mathcal{P}_h u_0 $$
is the spatial discretization error, while

$$e_2^n := (e^{-t_n L_h} - (I + \Delta t L_h)^{-n}) \mathcal{P}_h u_0$$

is the temporal discretization error. Clearly, we have $e_1^n = 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

$$\frac{\partial y(t)}{\partial t} + Ly(t) = 0, \quad y(0) = u_0$$

and its finite element semi-discrete equation

$$\frac{\partial y_h(t)}{\partial t} + L_h y_h(t) = 0, \quad y_h(0) = \mathcal{P}_h u_0$$

respectively. Let

$$e_1(t) := y(t) - y_h(t), \quad \rho(t) := L_h^{-1} \frac{\partial e_1(t)}{\partial t} + e_1(t).$$

It can be verified that

$$\rho(t) = (L^{-1} - L_h^{-1}) Ly(t).$$

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$,

$$\|e_1(t)\|_{L^2(D)} \le c \sup_{0 \le \tau \le t} \left( \|\rho(\tau)\|_{L^2(D)} + \tau \left\| \frac{\partial \rho(\tau)}{\partial \tau} \right\|_{L^2(D)} \right).$$

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:

$$\|\rho(\tau)\|_{L^2(D)} = \|(L^{-1} - L_h^{-1})Ly(\tau)\|_{L^2(D)} \le c h^2 \|Le^{-L\tau}u_0\|_{L^2(D)} \le c h^2,$$

$$\tau \left\| \frac{\partial \rho(\tau)}{\partial \tau} \right\|_{L^2(D)} = \tau \left\| (L^{-1} - L_h^{-1}) L^2 y(\tau) \right\|_{L^2(D)} \le c \tau h^2 \left\| L^2 e^{-L\tau} u_0 \right\|_{L^2(D)} \le c h^2 \|Lu_0\|_{L^2(D)} \le c h^2.$$

Thus

$$\|e_1^n\|_{L^2(D)} = \|e_1(t_n)\|_{L^2(D)} \le c h^2, \quad n = 0, 1, \dots, N. \qquad (4.13)$$

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

$$\begin{aligned} \|e_2^n\|_{L^2(D)} &= \left\| (e^{-n\Delta t L_h} - (I + \Delta t L_h)^{-n}) L_h^{-1} L_h \mathcal{P}_h u_0 \right\|_{L^2(D)} \\ &\le \left\| (e^{-n\Delta t L_h} - (I + \Delta t L_h)^{-n}) L_h^{-1} \right\|_{\mathcal{L}(L^2(D))} \|L_h P_h u_0\|_{L^2(D)} \end{aligned}$$

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

$$\sup_{j=1,...,K} \left| \frac{e^{-n\Delta t \lambda_j^h} - (1 + \Delta t \lambda_j^h)^{-n}}{\lambda_j^h} \right|,$$
where $\lambda_j^h > 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

$$ \sup_{j=1,...,K} \left| \frac{(e^{-n\Delta t}\lambda_j^h - (1+\Delta t\lambda_j^h)^{-n})}{\lambda_j^h} \right| \le c\Delta t. $$

This proves

$$ \|e_2^n\|_{L^2(D)} \le c\Delta t, \quad n = 0, 1, \dots, N. \qquad (4.14) $$

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

$$ \|\theta_1\|_{L^2(D)} \le c(\Delta t + h^2). $$

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

$$ \|\theta_1\|_{L^2(\Omega, L^2(D))} \le c(\Delta t + h^2). $$

**Remark 4.1.** If the $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$,

$$ \|\theta_1\|_{L^2(D)} = \|T_n u_0\|_{L^2(D)} \le c\|u_0\|_{L^2(D)} \frac{\Delta t + h^2}{t_n}, \quad n = 1, \dots, N. \qquad (4.15) $$

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

$$ \|\theta_2\|_{L^2(\Omega, L^2(D))} \le c[\Delta t^{\frac{1}{2}} + (\Delta t + h^2)\ln(\Delta t^{-1}) + \sum_{k=0}^{n-1} \|u(t_k) - u_h^k\|_{L^2(\Omega, L^2(D))} \Delta t], \quad n = 1, \dots, N, $$

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

*Proof.* The term to be bounded can be expressed by

$$ \theta_2 = \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} (S(t_n - \tau)f(u(\tau)) - S_{h,\Delta t}^{n-k}\mathcal{P}_h f(u_h^k))d\tau, $$

which can be decomposed into

$$ \theta_2 = \theta_2^1 + \theta_2^2 + \theta_2^3 + \theta_2^4 $$

with

$$ \theta_2^1 = \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} (S(t_n - \tau) - S(t_n - t_k))f(u(\tau))d\tau, \quad \theta_2^2 = \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} (S(t_n - t_k) - S_{h,\Delta t}^{n-k}\mathcal{P}_h)f(u(\tau))d\tau, $$
$$ \theta_2^3 = \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} S_{h, \Delta t}^{n-k} \mathcal{P}_h(f(u(\tau)) - f(u(t_k))) d\tau, \quad \theta_2^4 = \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} S_{h, \Delta t}^{n-k} \mathcal{P}_h(f(u(t_k)) - f(u_h^k)) d\tau. $$

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

$$ \|\theta_2^1\|_{L^2(\Omega, L^2(D))} \le \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \mathbb{E} \left[ \|S(t_n - \tau) - S(t_n - t_k)\|_{L^2(D)}^2 \|f(u(\tau))\|_{L^2(D)}^2 \right]^{\frac{1}{2}} d\tau. $$

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

$$ \begin{aligned} \|S(t_n - \tau) - S(t_n - t_{n-1})\|_{L^2(D)} &\le c, \\ \|S(t_n - \tau) - S(t_n - t_k)\|_{L^2(D)} &\le \|LS(t_n - \tau)\|_{L^2(D)} \|L^{-1}(I - S(\tau - t_k))\|_{L^2(D)} \\ &\le c \frac{\tau - t_k}{t_n - \tau}, \quad \tau \in (t_k, t_{k+1}), \quad k = 0, \dots, n-2. \end{aligned} $$

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

$$ \begin{aligned} \|\theta_2^1\|_{L^2(\Omega, L^2(D))} &\le c(\Delta t + \sum_{k=0}^{n-2} \int_{t_k}^{t_{k+1}} \frac{\Delta t}{t_n - t_{k+1}} d\tau) \le c\Delta t(1 + \sum_{k=1}^{n} \frac{1}{k}) \\ &\le c\Delta t(1 + \ln n) \le c\Delta t(1 + \ln(\Delta t^{-1})). \end{aligned} $$

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

$$ \begin{aligned} \|\theta_2^2\|_{L^2(\Omega, L^2(D))} &\le \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \|T_{n-k}f(u(\tau))\|_{L^2(\Omega, L^2(D))} d\tau \\ &\le \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \mathbb{E}[\|T_{n-k}\|_{L^2(D)}^2 \|f(u(\tau))\|_{L^2(D)}^2]^{\frac{1}{2}} d\tau \\ &\le c \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \frac{\Delta t + h^2}{t_{n-k}} d\tau = c(\Delta t + h^2) \sum_{k=0}^{n-1} \frac{1}{n-k} \le c(\Delta t + h^2) \ln(\Delta t^{-1}). \end{aligned} $$

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

$$ \begin{aligned} \|\theta_2^3\|_{L^2(\Omega, L^2(D))} &\le \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \|f(u(\tau)) - f(u(t_k))\|_{L^2(\Omega, L^2(D))} d\tau \\ &\le c \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \|u(\tau) - u(t_k)\|_{L^2(\Omega, L^2(D))} d\tau \\ &\le c \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} (\tau - t_k)^{\frac{1}{2}} d\tau \le c\Delta t^{\frac{1}{2}}. \end{aligned} $$
The part $\theta_2^4$ can be estimated similarly:

$$
\begin{align*}
\|\theta_2^4\|_{L^2(\Omega, L^2(D))} &\le \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \|f(u(t_k)) - f(u_h^k)\|_{L^2(\Omega, L^2(D))} d\tau \\
&\le c \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \|u(t_k) - u_h^k\|_{L^2(\Omega, L^2(D))} d\tau \\
&= c \sum_{k=0}^{n-1} \|u(t_k) - u_h^k\|_{L^2(\Omega, L^2(D))} \Delta t.
\end{align*}
$$

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 = O(j^{-(2\gamma+1+\epsilon)})$ for some $\gamma \ge 0$ and $\epsilon > 0$. Then it holds: for $0 \le \tau_1 \le \tau_2 \le T$,

$$
\| \mathcal{P}_h (G(u(\tau_2)) - G(u(\tau_1)) \mathcal{P}_J^w ) \|_{L^2(\Omega, \mathcal{L}_Q)} \le c (|\tau_2 - \tau_1|^{\frac{1}{2}} + J^{-\gamma}). \quad (4.16)
$$

*Proof.* Using the triangle inequality:

$$
\| \mathcal{P}_h (G(u(\tau_2)) - G(u(\tau_1)) \mathcal{P}_J^w ) \|_{L^2(\Omega, \mathcal{L}_Q)} \le \| \mathcal{P}_h (G(u(\tau_2)) - G(u(\tau_1))) \|_{L^2(\Omega, \mathcal{L}_Q)} \\
+ \| \mathcal{P}_h (G(u(\tau_1)) - G(u(\tau_1)) \mathcal{P}_J^w ) \|_{L^2(\Omega, \mathcal{L}_Q)},
$$

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

$$
\begin{align*}
\| \mathcal{P}_h (G(u(\tau_2)) - G(u(\tau_1))) \|_{L^2(\Omega, \mathcal{L}_Q)} &\le \| G(u(\tau_2)) - G(u(\tau_1)) \|_{L^2(\Omega, \mathcal{L}_Q)} \\
&\le c \| u(\tau_2) - u(\tau_1) \|_{L^2(\Omega, L^2(D))} \le c |\tau_2 - \tau_1|^{\frac{1}{2}}.
\end{align*}
$$

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

$$
\begin{align*}
\| \mathcal{P}_h (G(u(\tau_1)) - G(u(\tau_1)) \mathcal{P}_J^w ) \|_{L^2(\Omega, \mathcal{L}_Q)} &\le \mathbb{E} \left[ \sum_{j=1}^{\infty} \| G(u(\tau_1))(I - \mathcal{P}_J^w) Q^{\frac{1}{2}} \phi_j \|_{L^2(D)}^2 \right]^{\frac{1}{2}} \\
&\le c \mathbb{E} \left[ \| G(u(\tau_1)) \|_{L^2(D)}^2 \sum_{j=1}^{\infty} \| (I - \mathcal{P}_J^w) q_j^{\frac{1}{2}} \phi_j \|_{L^2(D)}^2 \right]^{\frac{1}{2}} \\
&\le c (1 + \| u_0 \|_{L^2(\Omega, L^2(D))}) \left( \sum_{j=J+1}^{\infty} q_j \right)^{1/2} \\
&\le c J^{-\gamma}.
\end{align*}
$$

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

$$ \|\theta_3\|_{L^2(\Omega, L^2(D))}^2 \le c[(\Delta t^{\frac{1}{2}} + h^2)^2 + \sum_{k=0}^{n-1} \|u(t_k) - u_h^k\|_{L^2(\Omega, L^2(D))}^2 \Delta t], \quad n = 1, \dots, N, \quad (4.17) $$

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

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

$$ \theta_3^i := \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} X_i dW(\tau) $$

with

$$ \begin{aligned} X_1 &:= (S(t_n - \tau) - S(t_n - t_k))G(u(\tau)), \quad X_2 := (S(t_n - t_k) - S_{h, \Delta t}^{n-k} \mathcal{P}_h)G(u(\tau)), \\ X_3 &:= S_{h, \Delta t}^{n-k} \mathcal{P}_h (G(u(\tau)) - G(u(t_k))) \mathcal{P}_j^w, \quad X_4 := S_{h, \Delta t}^{n-k} \mathcal{P}_h (G(u(t_k))) \mathcal{P}_j^w - G(u_h^k) \mathcal{P}_j^w. \end{aligned} $$

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

$$ \begin{aligned} \|\theta_3^1\|_{L^2(\Omega, L^2(D))}^2 &= \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \mathbb{E}[\|X_1\|_{\mathcal{L}_Q}^2] d\tau \\ &\le \int_{t_{n-1}}^{t_n} \mathbb{E}[\|S(t_n - \tau) - S(t_n - t_{n-1})\|_{L(L^2(D))}^2 \|G(u(\tau))\|_{\mathcal{L}_Q}^2] d\tau \\ &\quad + \sum_{k=0}^{n-2} \int_{t_k}^{t_{k+1}} \mathbb{E}[\|(S(t_n - \tau) - S(t_n - t_k))L^{-\frac{1}{2}}\|_{L(L^2(D))}^2 \|L^{\frac{1}{2}}G(u(\tau))\|_{\mathcal{L}_Q}^2] d\tau. \end{aligned} $$

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

$$ \int_{t_{n-1}}^{t_n} \mathbb{E}[\|S(t_n - \tau) - S(t_n - t_{n-1})\|_{L(L^2(D))}^2 \|G(u(\tau))\|_{\mathcal{L}_Q}^2] d\tau \le c \Delta t. $$

Then employing (4.8) and (4.9) gives

$$ \begin{aligned} \|(S(t_n - \tau) - S(t_n - t_k))L^{-\frac{1}{2}}\|_{L(L^2(D))}^2 &= \|L^{\frac{1}{2}}S(t_n - \tau)L^{-1}(I - S(\tau - t_k))\|_{L(L^2(D))}^2 \\ &\le \|L^{\frac{1}{2}}S(t_n - \tau)\|_{L(L^2(D))}^2 \|L^{-1}(I - S(\tau - t_k))\|_{L(L^2(D))}^2 \\ &\le c \frac{(\tau - t_k)^2}{t_n - \tau}. \end{aligned} $$

Making use of (2.8), (2.11) gives

$$ \begin{aligned} & \sum_{k=0}^{n-2} \int_{t_k}^{t_{k+1}} \mathbb{E}[\|(S(t_n - \tau) - S(t_n - t_k))L^{-\frac{1}{2}}\|_{L(L^2(D))}^2 \|L^{\frac{1}{2}}G(u(\tau))\|_{\mathcal{L}_Q}^2] d\tau \\ &\le c \sum_{k=0}^{n-2} \int_{t_k}^{t_{k+1}} \frac{\Delta t^2}{t_n - t_{k+1}} d\tau \le c \Delta t^2 \ln(\Delta t^{-1}). \end{aligned} $$
Therefore

$$ \|\theta_3^1\|_{L^2(\Omega, L^2(D))}^2 \le c(\Delta t + \Delta t^2 \ln(\Delta t^{-1})). $$

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

$$ \begin{aligned} \|\theta_3^2\|_{L^2(\Omega, L^2(D))}^2 &= \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \mathbb{E}[\|X_2\|_{L_Q}^2] d\tau = \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \mathbb{E}[\|T_{n-k}G(u(\tau))\|_{L_Q}^2] d\tau \\ &\le c \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \mathbb{E}\left[(\frac{\Delta t + h^2}{t_{n-k}})^2 \|G(u(\tau))\|_{L_Q}^2\right] d\tau. \end{aligned} $$

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

$$ \|\theta_3^2\|_{L^2(\Omega, L^2(D))}^2 \le c \frac{(\Delta t + h^2)^2}{\Delta t} \sum_{k=0}^{n-1} \frac{1}{(n-k)^2} \le c \frac{(\Delta t + h^2)^2}{\Delta t} \le c\Delta t. $$

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

$$ \begin{aligned} \|\theta_3^3\|_{L^2(\Omega, L^2(D))}^2 &= \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \mathbb{E}[\|X_3\|_{L_Q}^2] d\tau \le c \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} (|\tau - t_k|^{\frac{1}{2}} + J^{-\gamma})^2 d\tau \\ &\le c \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} (|\tau - t_k|^{\frac{1}{2}} + h^2)^2 d\tau \\ &\le c(\Delta t^{\frac{1}{2}} + h^2)^2. \end{aligned} $$

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

$$ \begin{aligned} \|\theta_3^4\|_{L^2(\Omega, L^2(D))}^2 &= \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \mathbb{E}[\|X_4\|_{L_Q}^2] d\tau = \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \mathbb{E}[\|S_{h,\Delta t}^{n-k}\mathcal{P}_h(G(u(t_k))\mathcal{P}_j^w - G(u_h^k)\mathcal{P}_j^w)\|_{L_Q}^2] d\tau \\ &\le c \sum_{k=0}^{n-1} \int_{t_k}^{t_{k+1}} \mathbb{E}[\|u(t_k) - u_h^k\|_{L^2(D)}^2] d\tau \\ &= c \sum_{k=0}^{n-1} \|u(t_k) - u_h^k\|_{L^2(\Omega, L^2(D))}^2 \Delta t. \end{aligned} $$

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

$$ \|u(t_n) - u_h^n\|_{L^2(\Omega, L^2(D))} \le c(\Delta t^{\frac{1}{2}} + (\Delta t + h^2) \ln(\Delta t^{-1})), \quad n = 1, \dots, N. $$
**Proof.** It follows from (4.3), Lemmas 4.1-4.4, and the triangle inequality:

$$ \|u(t_n) - u_h^n\|_{L^2(\Omega, L^2(D))}^2 \le c\left[(\Delta t^{\frac{1}{2}} + (\Delta t + h^2) \ln(\Delta t^{-1}))^2 + \sum_{k=0}^{n-1} \|u(t_k) - u_h^k\|_{L^2(\Omega, L^2(D))}^2 \Delta t\right], n = 1, \dots, N. $$

Then the discrete Gronwall inequality yields

$$ \|u(t_n) - u_h^n\|_{L^2(\Omega, L^2(D))}^2 \le c(\Delta t^{\frac{1}{2}} + (\Delta t + h^2) \ln(\Delta t^{-1}))^2, \quad n = 1, \dots, N. $$

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

$$ \|u(t_n) - u_h^n\|_{L^2(\Omega, L^2(D))} \le c(\Delta t^{\frac{1}{2}} + h^2 \ln(\Delta t^{-1})), \quad n = 1, \dots, N. $$

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

$$ \|u(t_n) - u_h^n\|_{L^2(\Omega, L^2(D))} \le c(\Delta t^{\frac{1}{2}} + h^2 \Delta t^{-\epsilon_0}), \quad n = 1, \dots, N, $$

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

$$ \|u(t_n) - u_h^n\|_{L^2(\Omega, L^2(D))} \le c(\Delta t^{\frac{1}{2}} + h^{2-\epsilon_0}), \quad n = 1, \dots, N. $$

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

$$ \begin{aligned} du(x,t) &= \varepsilon \partial_x (e^{z(x,\omega)} \partial_x u) dt - \frac{u}{1+u^2} dt + G(u) dW(x,t), \quad 0 < t < T, \quad x \in (0,1), \\ u(0,t) &= u(1,t) = 0, \quad 0 \le t \le T, \\ u(x,0) &= u_0(x), \quad x \in [0,1], \end{aligned} \tag{5.1} $$

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_0^\gamma$-valued Wiener process defined by

$$ W(t, x) = \sum_{j=1}^{\infty} \sqrt{q_j} \sin(j\pi x) \beta_j(t), $$

where $q_j = 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 $z = (z(x_1, \omega), \dots, z(x_p, \omega))^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 $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 $z$ will gradually approach the theoretical mean-zero function. For each $x \in \bar{D}$, we denote by $\mu_m(x)$ the approximation of $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^{-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_j^{\text{ref}}$ 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_h^N\|_{L^2(\Omega, L^2(D))}$ is approximately calculated by

$$ \|u(T) - u_h^N\|_{L^2(\Omega, L^2(D))} \approx \left( \frac{1}{200} \sum_{j=1}^{200} \|u_j^{\text{ref}} - u_{j,h}^N\|^2_{L^2(D)} \right)^{\frac{1}{2}} =: u_{\text{error}}. $$

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_{\text{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_{\text{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.

<table><caption>**Table 1.** $\rho(\Lambda_{-})$ changes with respect to $M$.</caption><thead><tr><th><i>M</i></th><th>100</th><th>500</th><th>1000</th><th>2000</th><th>3000</th><th>4000</th></tr></thead><tbody><tr><th><i>ρ</i>(Λ<sub>-</sub>)</th><td>5.19</td><td>7.26E-01</td><td>4.27E-03</td><td>2.68E-08</td><td>1.58E-13</td><td>3.14E-14</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Mean of sample values $|\mu_m(x)|$ as function of $m$ in log-log scale.</figcaption></figure>

<table><caption>**Table 2.** Time accuracy test.</caption><thead><tr><th>Δt</th><th>u<sub>error</sub></th><th>Order</th></tr></thead><tbody><tr><td>1.00E-2</td><td>4.00E-3</td><td>—</td></tr><tr><td>5.00E-3</td><td>2.73E-3</td><td>0.55</td></tr><tr><td>2.50E-3</td><td>1.86E-3</td><td>0.55</td></tr><tr><td>1.25E-3</td><td>1.30E-3</td><td>0.52</td></tr><tr><td>6.25E-4</td><td>9.21E-4</td><td>0.50</td></tr></tbody></table>

<table><caption>**Table 3.** Time accuracy test.</caption><thead><tr><th>Δt</th><th>u<sub>error</sub></th><th>Order</th></tr></thead><tbody><tr><td>1.00E-2</td><td>5.17E-3</td><td>—</td></tr><tr><td>5.00E-3</td><td>3.81E-3</td><td>0.44</td></tr><tr><td>2.50E-3</td><td>2.82E-3</td><td>0.44</td></tr><tr><td>1.25E-3</td><td>1.90E-3</td><td>0.57</td></tr><tr><td>6.25E-4</td><td>1.31E-3</td><td>0.54</td></tr></tbody></table>

<table><caption>**Table 4.** Spatial accuracy test.</caption><thead><tr><th>h</th><th>u<sub>error</sub></th><th>Order</th></tr></thead><tbody><tr><td>1/16</td><td>1.20E-2</td><td>—</td></tr><tr><td>1/32</td><td>3.38E-3</td><td>1.83</td></tr><tr><td>1/64</td><td>9.76E-4</td><td>1.79</td></tr><tr><td>1/128</td><td>2.84E-4</td><td>1.78</td></tr><tr><td>1/256</td><td>7.32E-5</td><td>1.95</td></tr></tbody></table>

<table><caption>**Table 5.** Spatial accuracy test.</caption><thead><tr><th>h</th><th>u<sub>error</sub></th><th>Order</th></tr></thead><tbody><tr><td>1/16</td><td>1.15E-2</td><td>—</td></tr><tr><td>1/32</td><td>3.17E-3</td><td>1.85</td></tr><tr><td>1/64</td><td>9.30E-4</td><td>1.77</td></tr><tr><td>1/128</td><td>2.41E-4</td><td>1.95</td></tr><tr><td>1/256</td><td>6.48E-5</td><td>1.89</td></tr></tbody></table>

**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:

$$u_t(x, t) = \varepsilon \partial_{xx} u + u - u^3, \quad 0 < t < T, \ x \in (0, 1),$$

$$u(0, t) = u(1, t) = 0, \quad 0 \le t \le T,$$

$$u(x, 0) = u_0(x), \ x \in [0, 1].$$

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).

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Time evolution of numerical solution with different diffusion coefficients. (a): deterministic case, (b): random case with $q = 0.1$, (c): random case with $q = 2$.</figcaption></figure>

**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:

$$ 
\begin{aligned}
u_t(x, t) &= \Delta u + \frac{u - u^3}{\varepsilon^2}, && t \in (0, T), x \in (0, 1)^2, \\
du(x, t) &= \nabla \cdot (e^{z(x, \omega)} \nabla u) dt + \frac{u - u^3}{\varepsilon^2} dt + G(u) dW(x, t),
\end{aligned} 
\quad (5.2) $$

where $z(x, \omega)$ stands for a bounded two-dimensional Gaussian random field with mean-zero and
<figure><img src="image_6.png" /><figcaption>Figure 3. Time evolution of numerical solution with different noise. (a): deterministic case, (b): random case with $\gamma = 0.5$, (c): random case with $\gamma = 1$.</figcaption></figure>

covariance function $c_q(x)$. Let $x = (x_1, x_2)^T$, $x_1, x_2 \in (0, 1)$. Denote $W(x, t)$ by

$$W(x, t) := W(x_1, x_2, t) = \sum_{j_1, j_2=1}^{\infty} \sqrt{q_{j_1, j_2}} \sin(j_1 \pi x_1) \sin(j_2 \pi x_2) \beta_{j_1, j_2}(t),$$

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 \to 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:

$$V = -\frac{1}{R} = -\kappa, \qquad (5.3)$$

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 - 2t}$ due to $V(t) = \frac{dR(t)}{dt} = -\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×512 mesh and $\varepsilon = 7.5 \times 10^{-4}$, $\Delta t = 5 \times 10^{-5}$. The initial configuration is defined by:

$$\begin{cases} u(x, 0) = u(x_1, x_2, 0) = \tanh \frac{1.5+1.2 \cos(6\theta)-2\pi r}{\sqrt{2\varepsilon}}, \\ \theta = \arctan \frac{x_2-0.5}{x_1-0.5}, \\ r = \sqrt{(x_1-0.5)^2 + (x_2-0.5)^2}. \end{cases}$$

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.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Evolution of a star-shaped interface for $t = 0, 0.001, 0.005$. The first row: deterministic case. The second row: random case with $q = 0.5$ and $G(u) = 0$. The third row: random case with $q = 0.5$ and $G(u) = 10(1 - |u|)$.</figcaption></figure>

**6. Conclusions**

In this paper we study a stochastic evolution equation driven by a bounded log-Whittle-Matern 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.

References

1. E. J. Allen, S. J. Novosel, Z. Zhang, Finite element and difference approximation of some linear stochastic partial differential equations, Stochastics, 64 (1998), 117–142.
https://doi.org/10.1080/17442509808834159

2. A. Andersson, S. Larsson, Weak convergence for a spatial approximation of the nonlinear stochastic heat equation, Math. Comput., 85 (2016), 1335–1358.
https://doi.org/10.1090/mcom/3016

3. I. Babuška, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal., 45 (2007), 1005–1034.
https://doi.org/10.1137/050645142

4. I. Babuška, R. Tempone, G. E. Zouraris, Galerkin finite element approximations of stochastic elliptic partial differential equations, SIAM J. Numer. Anal., 42 (2004), 800–825.
https://doi.org/10.1137/S0036142902418680

5. M. Beccari, M. Hutzenthaler, A. Jentzen, R. Kurniawan, F. Lindner, D. Salimova, Strong and weak divergence of exponential and linear-implicit euler approximations for stochastic partial differential equations with superlinearly growing nonlinearities, arXiv preprint arXiv:1903.06066, 2019.

6. C. E. Bréhier, Approximation of the invariant measure with an Euler scheme for stochastic PDEs driven by space-time white noise, Potential Anal., 40 (2014), 1–40.

7. C. E. Bréhier, J. Cui, J. Hong, Strong convergence rates of semidiscrete splitting approximations for the stochastic Allen–Cahn equation, IMA J. Numer. Anal., 39 (2019), 2096–2134.
https://doi.org/10.1093/imanum/dry052

8. M. Cai, S. Gan, X. Wang, Weak convergence rates for an explicit full-discretization of stochastic Allen-Cahn equation with additive noise, J. Sci. Comput., 86 (2021), 1–30.

9. Y. Cao, J. Hong, Z. Liu, Approximating stochastic evolution equations with additive white and rough noises, SIAM J. Numer. Anal., 55 (2017), 1958–1981. https://doi.org/10.1137/16M1056122

10. N. Chopin, Fast simulation of truncated Gaussian distributions, Stat. Comput., 21 (2011), 275–288.
https://doi.org/10.1007/s11222-009-9168-1

11. J. Cui, J. Hong, Analysis of a splitting scheme for damped stochastic nonlinear Schrödinger equation with multiplicative noise, SIAM J. Numer. Anal., 56 (2018), 2045–2069.
https://doi.org/10.1137/17M1154904

12. J. Cui, J. Hong, Strong and weak convergence rates of a spatial approximation for stochastic partial differential equation with one-sided Lipschitz coefficient, SIAM J. Numer. Anal., 57 (2019), 1815–1841. https://doi.org/10.1137/18M1215554
13. J. Cui, J. Hong, Z. Liu, Strong convergence rate of finite difference approximations for stochastic cubic Schrödinger equations, *J. Differ. Equations*, **263** (2017), 3687–3713. https://doi.org/10.1016/j.jde.2017.05.002

14. J. Cui, J. Hong, Z. Liu, W. Zhou, Strong convergence rate of splitting schemes for stochastic nonlinear Schrödinger equations, *J. Differ. Equations*, **266** (2019), 5625–5663. https://doi.org/10.1016/j.jde.2018.10.034

15. J. Cui, J. Hong, L. Sun, Strong Convergence of Full Discretization for Stochastic Cahn-Hilliard Equation Driven by Additive Noise, *SIAM J. Numer. Anal.*, **59** (2021), 2866–2899. https://doi.org/10.1137/20M1382131

16. J. Cui, J. Hong, L. Sun, Weak convergence and invariant measure of a full discretization for parabolic SPDEs with non-globally Lipschitz coefficients, *Stoch. Proc. their Appl.*, **134** (2021), 55–93.

17. G. Da Prato, J. Zabczyk, Stochastic equations in infinite dimensions, *Cambridge university press*, 2014.

18. A. Davie, J. Gaines, Convergence of numerical schemes for the solution of parabolic stochastic partial differential equations, *Math. Comput.*, **70** (2001), 121–134. https://doi.org/10.1090/S0025-5718-00-01224-2

19. A. De Bouard, A. Debussche, Weak and strong order of convergence of a semidiscrete scheme for the stochastic nonlinear Schrödinger equation, *Appl. Math. Optim.*, **54** (2006), 369–399. https://doi.org/10.1007/s00245-006-0875-0

20. A. Debussche, Weak approximation of stochastic partial differential equations: The nonlinear case, *Math. Comput.*, **80** (2011), 89–117. https://doi.org/10.1090/S0025-5718-2010-02395-6

21. A. Debussche, J. Printems, Weak order for the discretization of the stochastic heat equation, *Math. Comput.*, **78** (2009), 845–863. https://doi.org/10.1090/S0025-5718-08-02184-4

22. C. R. Dietrich, A simple and efficient space domain implementation of the turning bands method, *Water Resour. Res.*, **31** (1995), 147–156. https://doi.org/10.1029/94WR01457

23. C. R. Dietrich, G. N. Newsam, Fast and exact simulation of stationary Gaussian processes through circulant embedding of the covariance matrix, *SIAM J. Sci. Comput.*, **18** (1997), 1088–1107. https://doi.org/10.1137/S1064827592240555

24. Q. Du, T. Zhang, Numerical approximation of some linear stochastic partial differential equations driven by special additive noises, *SIAM J. Numer. Anal.*, **140** (2002), 1421–1445.

25. K. Engel, R. Nagel, One-parameter semigroups for linear evolution equations, *Semigroup Forum*, **63** (1999), 278–280.

26. X. Feng, Y. Li, Y. Zhang, Finite element methods for the stochastic Allen-Cahn equation with gradient-type multiplicative noise, *SIAM J. Numer. Anal.*, **55** (2017), 194–216. https://doi.org/10.1137/15M1022124

27. M. Geissert, M. Kovács, S. Larsson, Rate of weak convergence of the finite element method for the stochastic heat equation with additive noise, *BIT Numer. Math.*, **549** (2009), 343–356.

28. I. Gohberg, S. Goldberg, M. A. Kaashoek, Hilbert-schmidt operators, In: *Classes of Linear Operators Vol. I*, 138–147. Springer, 1990.
29. I. Gyöngy, S. Sabanis, D. Šiška, Convergence of tamed Euler schemes for a class of stochastic evolution equations, *Stochast. Partial Differ. Equations: Anal. Comput.*, **4** (2016), 225–245.

30. I. Gyöngy, A. Millet, Rate of convergence of space time approximations for stochastic evolution equations, *Potential Anal.*, **30** (2009), 29–64. https://doi.org/10.1111/j.1755-3768.1986.tb06988.x

31. E. Hausenblas, Approximation for semilinear stochastic evolution equations, *Potential Anal.*, **18** (2003), 141–186. https://doi.org/10.1023/A:1020552804087

32. E. Hausenblas, Weak approximation for semilinear stochastic evolution equations, In: *Stochastic analysis and related topics VIII*, 111–128. Springer, 2003.

33. M. Hutzenthaler, A. Jentzen, *Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients*, volume 236. American Mathematical Society, 2015.

34. M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients, *P. Royal Soc. A-Math. Phy.*, **467** (2011), 1563–1576. https://doi.org/10.1098/rspa.2010.0348

35. A. Jentzen, P. E. Kloeden, Overcoming the order barrier in the numerical approximation of stochastic partial differential equations with additive space-time noise, *P. Royal Soc. A-Math. Phy.*, **465** (2008), 649–667. https://doi.org/10.1098/rspa.2008.0325

36. A. Jentzen, P. Pušnik, Strong convergence rates for an explicit numerical approximation method for stochastic evolution equations with non-globally Lipschitz continuous nonlinearities, *IMA J. Numer. Anal.*, **40** (2020), 1005–1050.

37. N. L. Johnson, S. Kotz, N. Balakrishnan, *Continuous univariate distributions*, volume 289, John wiley & sons, 1995.

38. Y. Kazashi, Quasi-monte carlo integration with product weights for elliptic PDEs with log-normal coefficients, *IMA J. Numer. Anal.*, **39** (2019), 1563–1593. https://doi.org/10.1093/imanum/dry028

39. P. E. Kloeden, G. J. Lord, A. et al Neuenkirch, The exponential integrator scheme for stochastic partial differential equations: Pathwise error bounds, *J. Comput. Appl. Math.*, **235** (2011), 1245–1260. https://doi.org/10.1016/j.cam.2010.08.011

40. P. E. Kloeden, E. Platen, Numerical solution of stochastic differential equations, *Springer Science and Business Media*, 2013.

41. M. Kovács, S. Larsson, F. Lindgren, Weak convergence of finite element approximations of linear stochastic evolution equations with additive noise, *BIT Numer. Math.*, **52** (2012), 85–108. https://doi.org/10.1007/s10543-011-0344-2

42. M. Kovács, S. Larsson, F. Lindgren, Weak convergence of finite element approximations of linear stochastic evolution equations with additive noise II, Fully discrete schemes, *BIT Numer. Math.*, **53** (2013), 497–525.

43. M. Kovács, S. Larsson, F. Lindgren, On the discretisation in time of the stochastic Allen-Cahn equation, *Math. Nachr.*, **291** (2018), 966–995. https://doi.org/10.1002/mana.201600283

44. R. Kruse, Optimal error estimates of Galerkin finite element methods for stochastic partial differential equations with multiplicative noise, *IMA J. Numer. Anal.*, **34** (2014), 217–251. https://doi.org/10.1093/imanum/drs055
45. R. Kruse, *Strong and weak approximation of semilinear stochastic evolution equations*, Springer, 2014.

46. C. Li, Y. Huang, N. Yi, An unconditionally energy stable second order finite element method for solving the Allen–Cahn equation, *J. Comput. Appl. Math.*, **353** (2019), 38–48.

47. N. Li, B. Meng, X. Feng, D. Gui, The spectral collocation method for the stochastic Allen-Cahn equation via generalized polynomial chaos, *Numer. Heat Tr. B-Fund.*, **68** (2015), 11–29.

48. N. Li, J. Zhao, X. Feng, D. Gui, Generalized polynomial chaos for the convection–diffusion equation with uncertainty, *Int. J. Heat Mass Transfer*, **97** (2016), 289–300. https://doi.org/10.1016/j.ijheatmasstransfer.2016.02.006

49. Y. Li, H. G. Lee, D. Jeong, J. Kim, An unconditionally stable hybrid numerical method for solving the Allen–Cahn equation, *Comput. Math. Appl.*, **60** (2010), 1591–1606. https://doi.org/10.1016/j.camwa.2010.06.041

50. F. Lindner, R. Schilling, Weak order for the discretization of the stochastic heat equation driven by impulsive noise, *Potential Anal.*, **38** (2013), 345–379. https://doi.org/10.1007/s11118-012-9276-y

51. Z. Liu, Z. Qiao, Strong approximation of monotone stochastic partial differential equations driven by white noise, *IMA J. Numer. Anal.*, **40** (2020), 1074–1093. https://doi.org/10.1093/imanum/dry088

52. Z. Liu, Z. Qiao, Strong approximation of monotone stochastic partial differential equations driven by multiplicative noise, *Stoch. Partial Differ.*, **9** (2021), 559–602. https://doi.org/10.1007/s40072-020-00179-2

53. G. J. Lord, C. E. Powell, T. Shardlow, An introduction to computational stochastic PDEs, *Cambridge University Press*, 2014.

54. A. Mantoglou, J. L. Wilson, The turning bands method for simulation of random fields using line generation by a spectral method, *Water Resour. Res.*, **18** (1982), 1379–1394. https://doi.org/10.1029/WR018i005p01379

55. G. Milstein, M. V. Tretyakov, Numerical integration of stochastic differential equations with nonglobally Lipschitz coefficients, *SIAM J. Numer. Anal.*, **43** (2005), 1139–1154. https://doi.org/10.1137/040612026

56. G. N. Milstein, M. V. Tretyakov, Stochastic numerics for mathematical physics, *Springer Science and Business Media*, 2013.

57. G. N. Newsam, C. R. Dietrich, Bounds on the size of nonnegative definite circulant embeddings of positive definite Toeplitz matrices, *IEEE T. Inform. Theory*, **40** (1994), 1218–1220. https://doi.org/10.1109/18.335952

58. J. Printems, On the discretization in time of parabolic stochastic partial differential equations, *ESAIM: Math. Model. Numer. Anal.*, **35** (2001), 1055–1078. https://doi.org/10.1051/m2an:2001148

59. M. Sauer, W. Stannat, Lattice approximation for stochastic reaction diffusion equations with one-sided Lipschitz condition, *Math. Comput.*, **84** (2015), 743–766. https://doi.org/10.1090/S0025-5718-2014-02873-1
60. M. Shinozuka, Simulation of multivariate and multidimensional random processes, *J. Acoust. Soc. Am.*, **49** (1971), 357–368. https://doi.org/10.1121/1.1912338

61. M. Shinozuka, C. M. Jan, Digital simulation of random processes and its applications, *J. Sound Vib.*, **25** (1972), 111–128. https://doi.org/10.1016/0022-460X(72)90600-1

62. J. L. Wadsworth, J. A. Tawn, Efficient inference for spatial extreme value processes associated to log-Gaussian random functions, *Biometrika*, **101** (2014), 1–15.

63. J. B. Walsh, Finite element methods for parabolic stochastic PDEs, *Potential Anal.*, **23** (2005), 1–43.

64. X. Wang, Strong convergence rates of the linear implicit Euler method for the finite element discretization of SPDEs with additive noise *IMA J. Numer. Anal.*, **37** (2017), 965–984.

65. X. Wang, An efficient explicit full-discrete scheme for strong approximation of stochastic Allen-Cahn equation, *Stoch. Proc. Appl.*, **130** (2020), 6271–6299. https://doi.org/10.1016/j.spa.2020.05.011

66. X. Wang, S. Gan, Weak convergence analysis of the linear implicit Euler method for semilinear stochastic partial differential equations with additive noise, *J. Math. Anal. Appl.*, **398** (2013), 151–169. https://doi.org/10.1016/j.jmaa.2012.08.038

67. A. TA. Wood, G. Chan, Simulation of stationary Gaussian processes in $[0, 1]^d$. *J. Comput. Graph. Stat.*, **3** (1994), 409–432. https://doi.org/10.1086/174579

68. D. Xiu, J. Shen, Efficient stochastic Galerkin methods for random diffusion equations, *J. Comput. Phys.*, **228** (2009), 266–281. https://doi.org/10.1016/j.jcp.2008.09.008

69. Y. Yan, Galerkin finite element methods for stochastic parabolic partial differential equations, *SIAM J. Numer. Anal.*, **43** (2005), 1363–1384. https://doi.org/10.1137/040605278

70. Z. Zhang, G. Karniadakis, *Numerical methods for stochastic partial differential equations with white noise*, Springer, 2017.

<figure><img src="image_8.png" /></figure>

AIMS Press

©2022 the Author(s), licensee AIMS Press. This
is an open access article distributed under the
terms of the Creative Commons Attribution License
(http://creativecommons.org/licenses/by/4.0)