A fully discrete HDG ensemble Monte Carlo algorithm for a heat equation under uncertainty

1.   Introduction

This paper concentrates on the numerical simulation of a heat problem involving random diffusion and Robin coefficients. The problem is to find $u(t, \mathbf{x}, \omega)$ and $p(t, \mathbf{x}, \omega)$ such that

Here, $D \subseteq \mathbb{R}^{d}$ ($d = 2, 3$) is a Lipschitz domain with boundary $\partial D = \partial D_{0} \cup \partial D_{1}$, where $\partial D_{0}$ and $\partial D_{1}$ are two non-overlapping parts. The vector $\mathbf{n}$ is the unit outward normal to $\partial D$, and $\Omega$ denotes the sample space.

Many problems in engineering and physics involve uncertainty. For example, in heat transfer dynamics, material properties, Robin boundary conditions (convective heat transfer coefficient), thermal conductivity (diffusion coefficient), etc., are affected by various forms of uncertainty. These uncertainties make the accurate prediction and analysis more challenging. For these problems, many numerical approaches have been devised ^[1,2]. In addition to stochastic finite element methods ^[3], stochastic collocation methods ^[4,5,6,7], and polynomial chaos methods ^[8], the Monte Carlo (MC) method is widely regarded as another very important approach (see ^[9,10,11,12]). The advantages of the MC method are that its convergence rate is independent of the dimensionality of the random parameters, and it is easy to implement. For the MC method, we first perform $M$ independent and identically distributed (i.i.d.) samples of random variable $\omega$. The $i$-th $(i = 1, 2, \cdots, M)$ realization $\omega_i$ satisfies

where we assume that $\rho_{i}(t, \mathbf{x}) = \rho(t, \mathbf{x}, \omega_i)$ and $\kappa_{i}(t, \mathbf{x}) = \kappa(t, \mathbf{x}, \omega_i)$ for $i = 1, 2, \cdots, M$, with similar expressions to the other variables. The average

of the solution for problem (1.2) is used as an approximation of the solution for problem (1.1).

Although the MC method is simple and easy to implement, its convergence speed is very slow. To improve this method, Jiang and Layton ^[13] proposed an ensemble approach for the random evolution equation. Since its proposal, this method has been widely promoted and applied ^{[12,14,15,16,17,18]}. In ^[16], a parabolic problem with a random diffusion coefficient is solved by the ensemble MC and finite element method. The error estimate is not optimal in space. To improve this, Yong et al. in ^[19] presented an optimal error estimate using the finite element method with ensemble MC. For the same problem as ^[16], Li and Luo used the ensemble MC and HDG method to approximate it and obtain the optimal error estimate about space. Similarly, for a parabolic problem with random diffusion and Robin coefficients, Yao et al. ^[18] used the ensemble MC and finite element method to approximate it and find the sub-optimal error estimate in space.

The discontinuous Galerkin (DG) method is an excellent approximation method for problems concerned with partial differential equations (PDEs). The DG method is particularly well-suited for problems with discontinuous coefficients. However, the main disadvantage of DG methods is the large number of degrees of freedom. The HDG method (see, e.g., ^[20]) approximates the solution's flux and trace by introducing numerical fluxes and numerical traces in a mixed formulation. Compared to DG methods, HDG methods significantly reduce the number of globally coupled degrees of freedom. To date, the HDG method has been utilized for a variety of problems, including convection-diffusion problems, flow problems, and hyperbolic equations (see, e.g., ^{[15,21,22,23]}).

In this work, we study the numerical approximate of a parabolic problem with random diffusion and Robin coefficients by the ensemble MC and HDG methods. By introducing the flux $\mathit{\boldsymbol{p}}$, the parabolic problem with random Robin coefficients and diffusion coefficients can be transformed into problem (1.1). After using MC sampling, the problem involves solving a large number of problems (1.2). Introduce two averages for the Robin coefficient and diffusion coefficient, respectively, and construct an ensemble format. Then in the calculation of each time step, a coefficient matrix can be shared, which reduces the computational complexity (the details can be seen in Eq (3.12)). The FEMC-HDG method has first-order temporal accuracy and optimal $L^2$ convergence in spatial space.

The paper is structured as follows: In Section 2, the necessary notations and preliminaries are provided. In Section 3, we introduce the fully discrete HDG ensemble scheme for problem (1.2), along with a comprehensive analysis of its stability and convergence. In Section 4, two numerical examples are provided to illustrate the effectiveness of the proposed method. Finally, concluding remarks are presented.

2.   Preliminaries

In this section, we introduce some notations that will be used throughout this work.

Let $(\cdot, \cdot)$ and $\|\cdot\|$ be the inner product and the $L^2(D)$ norm, respectively. We adopt the standard Sobolev space notation $W^{s, q}(D)$, with the corresponding norm $\|\cdot\|_{W^{s, q}(D)}$ and seminorm $|\cdot|_{W^{s, q}(D)}$, where $s \geq 0$ and $1 \leq q \leq \infty$. For convenience, we use $H^{s}(D) : = W^{s, 2}(D)$. Specifically, the norm $\|\cdot\|_{H^{s}(D)}$ and semi-norm $|\cdot|_{H^{s}(D)}$ correspond to $H^{s}(D)$.

We assume that $(\Omega, \mathcal{F}, \mathbb{P})$ is a complete probability space, where $\mathcal{F} \subset 2^{\Omega}$ represents the $\sigma$-algebra of events, and $\mathbb{P}: \mathcal{F} \rightarrow [0, 1]$ denotes the probability measure. $H \in L_{\mathbb{P}}^{1}(\Omega)$ is a random variable. The expectation of $H$ can be expressed as

Consider a $d$-tuple $\delta = (\delta_{1}, \cdots, \delta_{d})$ with length $|\delta| = \sum_{i = 1}^{d} \delta_{i}$, where each $\delta_{i} \in \mathbb{N}^{+}$. Let the space $\widetilde{W}^{s, p}(D) = L_{\mathbb{P}}^{p}(\Omega, W^{s, p}(D))$ consist of random functions $u: \Omega \times D \to \mathbb{R}$, which are measurable with respect to the $\sigma$-algebra $\mathcal{F} \otimes \mathcal{B}(D)$, where $\mathcal{B}(D)$ denotes the Borel $\sigma$-algebra on $D$. The norm in $\widetilde{W}^{s, p}(D)$ is defined as

If $u \in \widetilde{W}^{s, p}(D)$, then for almost every $\omega$, $u(\omega, \cdot) \in W^{s, p}(D)$. Additionally, for $|\delta| \leq s$, the derivatives $\partial^{\delta} u(\cdot, \mathbf{x})$ are in $L_{\mathbb{P}}^{p}(\Omega)$ for almost every point on the domain $D$. Lastly, we define $\widetilde{H}^{s}(D) = L_{\mathbb{P}}^{2}(\Omega, H^{s}(D))$.

We introduce the tensor product Hilbert space

where the inner product is defined as

The norm is expressed as

Let $\mathcal{T}_{h}$ be a quasi-uniform triangulation of $D$. $\bar{\partial D_0}\cap\bar{\partial D_0}$ are the triangulation nodes. Thus, $D$ can be written as the union ${D} = \bigcup_{K \in \mathcal{T}_{h}} {K}$, where each element $K$ has a diameter $h_{K}$. We further define $h$ as $h = \max_{K \in \mathcal{T}_{h}} h_{K}$. The set $\partial \mathcal{T}_h$ is defined as $\{\partial K : K \in \mathcal{T}_h\}$, and $\mathcal{F}_h$ represents the collection of faces $F$ from elements $K \in \mathcal{T}_h$. The space of polynomials of degree $\ell$ on $K$ is denoted by $\mathbf{P}_{\ell}(K)$. We then introduce the following discontinuous finite element spaces:

The inner products are defined as follows:

where $(v, w)_D: = {\int}_D v wd\mathbf{x}$, and $\langle v, w \rangle_{\partial D}: = {\int}_{\partial D} v wds$. So we define $\|w\|_{\mathcal{T}_h}^2 = (w, w)_{\mathcal{T}_h}$.

We assume that the constant $C$ is positive and changes throughout the paper for different occurrences. Crucially, it does not depend on the discrete parameters $M$, $\Delta t$, or $h$.

3.   The FEMC-HDG scheme for the stochastic heat equation

In this section, we begin with outlining several assumptions, then develop an FEMC-HDG scheme for problem (1.2), and proceed to design the FEMC-HDG algorithm. Following this, we present results on stability and error estimates for both the FEMC-HDG scheme and the FEMC-HDG algorithms.

The mean values of the Robin boundary conditions and diffusion coefficients across the ensemble can be defined as follows:

and

respectively. In order to derive the stability and error estimates for the FEMC-HDG algorithm, we draw on the work of ^[16] and make the following two assumptions:

$(\mathbf{H1}):$ Let $\kappa_{\max}$, $\kappa_{\min}$, $\rho_{\min}$, and $\rho_{\max}$ be four positive constants such that for any $t \in [0, T]$, the probability is expressed as

and

$(\mathbf{H2}):$ Let $\kappa_{+}$, $\kappa_{-}$, $\rho_{+}$, and $\rho_{-}$ be four positive constants, such that for every $t \in [0, T]$, the probability is stated as follows:

and

Hypothesis ($\mathbf{H1}$) ensures almost sure uniform coercivity, while hypothesis ($\mathbf{H2}$) imposes a uniform bound on $|\kappa(t, \mathbf{x}, \omega_{i}) - \bar{\kappa}(t, \mathbf{x})|$ with a probability close to certainty. Similar properties are also assumed to hold for $\rho(t, \mathbf{x}, \omega_{i})$.

For a discretized physical space, problem (1.2) seeks $({\mathit{\boldsymbol{p}}}_{ih}(t, \cdot), u_{ih}(t, \cdot), \widehat u_{ih}(t, \cdot))\in{\mathit{\boldsymbol{V}}}^\ell_h \times W^\ell_h \times M^\ell_h$ such that, for all $({\mathit{\boldsymbol{r}}}_h, v_h, \mu_h) \in {\mathit{\boldsymbol{V}}}^\ell_h \times W^\ell_h \times M^\ell_h$,

The numerical flux $\widehat{{\mathit{\boldsymbol{p}}}}_{ih}$ is defined by:

where $\tau$ represents the stabilization parameter. In this work, we set $\tau = 1$ since we do not address multiple physical scales. Substituting Eq (3.10) into Eq (3.9), we obtain the semidiscrete HDG scheme of problem (1.2), which seeks $({\mathit{\boldsymbol{p}}}_{ih}(t, \cdot), u_{ih}(t, \cdot), \widehat u_{ih}(t, \cdot))\in {\mathit{\boldsymbol{V}}}^\ell_h\times W^\ell_h\times M^\ell_h$ such that, for all $({\mathit{\boldsymbol{r}}}_h, v_h, \mu_h) \in {\mathit{\boldsymbol{V}}}^\ell_h \times W^\ell_h \times M^\ell_h$, $i = 1, 2, \cdots, M,$

We divide the interval $[0, T]$ into $N$ equal parts, with a length of $\Delta t$ for each part. Denote $t_{n} = n \Delta t$ for $n = 1, 2, \cdots, N$. At each time $t = t_{n}$, the functions $u_{i}$, $f_{i}$, $g_{i}$, $\rho_{i}$, $\kappa_{i}$, $\bar{\rho}$, and $\bar{\kappa}$ are denoted by $u_{i}^{n}$, $f_{i}^{n}$, $g_{i}^{n}$, $\rho_{i}^{n}$, $\kappa_{i}^{n}$, $\bar{\rho}^{n}$, and $\bar{\kappa}^{n}$, respectively. Using the backward Euler method to approximate the time derivative of Eq (3.11), we can get the fully discrete non-ensemble MC HDG (FNEMC-HDG) scheme. Here, we omit the details. By applying Eqs (3.1) and (3.2) to (3.11) and performing simple algebraic calculations, we can obtain the FEMC-HDG scheme: Seek $({\mathit{\boldsymbol{p}}}^n_{ih}, u^n_{ih}, \widehat u^n_{ih})\in {\mathit{\boldsymbol{V}}}^\ell_h \times W^\ell_h \times M^\ell_h$ such that, for all $({\mathit{\boldsymbol{r}}}_h, v_h, \mu_h) \in {\mathit{\boldsymbol{V}}}^\ell_h \times W^\ell_h \times M^\ell_h$,

The initial conditions are given by $u_{ih}^{0} = \Pi^{\ell+1} u^0_i$, $p^0_{ih} = \dfrac{1}{\kappa_i} \nabla u_{ih}^{0}$, and $\widehat{u}^0_{ih} = \widehat{\Pi}^{\ell} u^0_i$, where $\Pi^{\ell+1}$ and $\widehat{\Pi}^{\ell}$ represent the $L^2$ projection operators $\Pi^{\ell+1}: L^2(K) \rightarrow \mathbf{P}_{\ell+1}(K)$ and $\widehat{\Pi}^{\ell}: L^2(F) \rightarrow \mathbf{P}_{\ell}(F)$, respectively, which satisfy:

To solve the stochastic partial differential equation (SPDE) (1.1) using the FEMC-HDG scheme, we first employ the MC method to generate i.i.d. samples. Subsequently, we apply the FEMC-HDG scheme to solve the resulting Eq (1.2). The solution of Eq (1.2) is used to compute the expected value of the solution to the SPDE (1.1). We present the FEMC-HDG algorithm which consists of three main steps:

Step 1. Generate an i.i.d. sample set for the initial conditions, boundary conditions, source term, diffusion coefficients, and Robin coefficients. For the $i$-th realization, these samples are denoted as $u_{i}^{0} = u^{0}\left(\cdot, \omega_{i}\right)$, $g_{i} = g\left(\cdot, \cdot, \omega_{i}\right)$, $f_{i} = f\left(\cdot, \cdot, \omega_{i}\right)$, $\kappa_{i} = \kappa\left(\cdot, \cdot, \omega_{i}\right)$, and $\rho_{i} = \rho\left(\cdot, \cdot, \omega_{i}\right)$, respectively.

Step 2. Calculate $\bar{\rho}^{n} = \frac{1}{M} \sum\limits_{i = 1}^{M} \rho\left(t_{n}, \mathbf{x}, \omega_{i}\right)$ and $\bar{\kappa}^{n} = \frac{1}{M} \sum\limits_{i = 1}^{M} \kappa\left(t_{n}, \mathbf{x}, \omega_{i}\right)$. For the $i$-th sample, solve Eq (3.12a) and (3.12b) to determine $u_{i, h}^{n}$ and ${\mathit{\boldsymbol{p}}}_{i, h}^{n}$ for $n = 1, \cdots, N$.

Step 3. For $n = 1, \cdots, N$, compute $\frac{1}{M} \sum\limits_{i = 1}^{M} u^n_{ih}$, $\frac{1}{M} \sum\limits_{i = 1}^{M} {\mathit{\boldsymbol{p}}}^n_{ih}$ to approximate the expectation $\mathbb{E}[u^n]$, $\mathbb{E}[{\mathit{\boldsymbol{p}}}^n]$, respectively.

The FEMC-HDG algorithm effectively combines the scheme (3.12a)–(3.12b) with the random sampling method. It preserves the benefits of the ensemble approach used for deterministic PDEs, such as employing the same coefficient matrix for all simulations at each time step. As a result, this approach involves solving a linear system with several right-hand side vectors, significantly reducing computational expenses (refer to ^[16]).

4.   Stability and convergence analysis

In this section, we carry out the stability analysis and error estimates.

4.1.   Stability

The FEMC-HDG scheme (3.12a)–(3.12b) has the following stability results.

Theorem 4.1. Assume that $f_{i} \in \widetilde{L}^{2}\left(0, T; L^2(D) \right)$, $g_{i} \in \widetilde{L}^{2}\left(0, T; L^{2}(\partial D_{1})\right)$, and $u^{0}_i\in \widetilde{L}^{2}\left(H^{1}(D)\right)$, and that hypotheses $(\mathbf{H1})$ and $(\mathbf{H2})$ are met. The FEMC-HDG scheme (3.12a)–(3.12b) is stable if

Moreover, the numerical solution to (3.12a)–(3.12b) satisfies the following inequality:

Proof. By selecting $(r_h, w_{h}, \mu_h) = (p^n_{ih}, u^n_{ih}, \widehat{u}_{ih}^n)$ in scheme (3.12a)–(3.12b), we obtain

Applying the polarization identity $(a-b) a = \frac{1}{2}[a^2 - b^2 + (a-b)^2]$ to Eqs (4.3) and (4.4), and then integrating over the probability space, we derive that

By utilizing the condition (3.6), and applying Young's inequality and the Cauchy-Schwarz inequality on the right-hand side of Eq (4.5), we obtain for any $\epsilon_i > 0$ $(i = 1, 2, 3)$ that

By combining the estimates given in Eqs (4.6) and (4.7), applying the conditions (3.3), (3.5), (3.7), and (3.8), and multiplying Eq (4.5) by $\Delta t$, we find that for all $0 < \sigma_1 \leq 1$,

According to the condition (4.1), choosing

from Eq (4.8), we obtain

Summing Eq (4.9) up from $n = 1$ to $n = N$, we get

Finally, applying Gronwall's inequality to the inequality in (4.10) yields the desired result.

Remark 4.1. According to condition (4.1), for the sequence $\left\{\kappa_{i}\right\}_{i = 1}^M$, the gap between $\kappa_{i}$ and $\bar{\kappa}$ should not exceed the fixed value of $\kappa_{min}$. A similar condition applies to the sequence $\left\{\rho_{i}\right\}_{i = 1}^M$. If this criterion is not satisfied, it becomes imperative to partition the ensemble into multiple smaller groups and employ the FEMC-HDG algorithm on each subset. Throughout this procedure, upholding the stability requirement within each subgroup is crucial.

4.2.   Error analysis

Under the assumption of sufficiently smooth solutions to the heat conduction equation, we provide error estimates for the FEMC-HDG algorithm. Let $\bar{U}_{Mh}^{n} = \frac{1}{M} \sum_{i = 1}^{M} u_{ih}^{n}$ and $\bar{{\mathit{\boldsymbol{Q}}}}_{Mh}^{n} = \frac{1}{M} \sum_{i = 1}^{M} {\mathit{\boldsymbol{p}}}_{ih}^{n}$ denote the results obtained from the FEMC-HDG algorithm, which approximate the expected exact solutions. We now proceed to estimate the error $\mathbb{E}\left[u_i{(t_{n})}\right]-\bar{U}_{Mh}^{n}$ and $\mathbb{E}\left[{\mathit{\boldsymbol{p}}}_i{(t_{n})}\right]-\bar{{\mathit{\boldsymbol{Q}}}}_{Mh}^{n}$ in various norms. For simplicity, we drop the subscript $i$ in $\mathbb{E}\left[u_i{(t_{n})}\right]-\bar{U}_{Mh}^{n}$ and $\mathbb{E}\left[{\mathit{\boldsymbol{p}}}_i{(t_{n})}\right]-\bar{{\mathit{\boldsymbol{Q}}}}_{Mh}^{n}$ and consider these errors in two separate terms:

The term $\mathcal{E}_{h}^{u^n} = \mathbb{E}\left[u(t_n)-u_{h}^{n}\right]$ represents the combined error due to HDG and time discretization. The second term, $\mathcal{E}_M^{u^n} = \mathbb{E}\left[u_{h}^{n}\right]-\bar{U}_{Mh}^{n}$, denotes the statistical error. $\mathcal{E}_{h}^{{\mathit{\boldsymbol{p}}}^n}$ and $\mathcal{E}^{{\mathit{\boldsymbol{p}}}^n}_{M}$ are analogous. Before deriving an error estimate for the FEMC-HDG algorithm, we undertake some preparatory work. First, we review standard error estimates for the $L^2$ projections ${\Pi}_{\ell}$ and $\widehat{\Pi}_{\ell}$.

Lemma 4.1. ^[24] Assume $\ell \geq 0$. Then there exists a constant $C$, which is not dependent on $K \in \mathcal{T}_h$, such that

Next, for any $t \in [0, T]$, let $(\mathit{\boldsymbol{\Pi}}^\ell_{V} \mathit{\boldsymbol{p}}_i, \Pi_{W}^\ell u_i)$ represent the HDG projection of $(\mathit{\boldsymbol{p}}_i, u_i)$. Here, $\mathit{\boldsymbol{\Pi}}_{V}^\ell \mathit{\boldsymbol{p}}_i$ and $\Pi_{W}^\ell u_i$ are the components of the HDG projection of $\mathit{\boldsymbol{p}}_i$ and $u_i$ into $\mathit{\boldsymbol{V}}^\ell_h$ and $W^\ell_h$, respectively. For each $K \in \mathcal{T}_h$, the pair $(\mathit{\boldsymbol{\Pi}}_{V}^\ell \mathit{\boldsymbol{p}}_i, \Pi_{W}^\ell u_i)$ satisfies the following equations:

for every face $F$ of the simplex $K$. When $\ell = 0$, the Eq (4.14a) and (4.14b) become trivial, and the projections are solely governed by Eq (4.14c). For these projections, they have the following property (Lemma 4.2), which is established in ^[20].

Lemma 4.2. Suppose the polynomial degree meets the requirements $\ell \geq 0$ and $\tau > 0$. Then the system (4.14) has a unique solution $(\mathbf{\Pi}^\ell_V {\mathit{\boldsymbol{p}}}_i, \Pi_W^\ell u_i)$. Moreover, there exists a constant $C > 0$, which is independent of both $K$ and $\tau$, such that

for $\ell_{u_i}$, $\ell_{{\mathit{\boldsymbol{p}}}_i}$ in $[0, \ell]$.

Lemma 4.3. For every $n = 1, 2, \cdots, N$, the subsequent equations:

and

hold for all $({\mathit{\boldsymbol{r}}}_h, v_h, \mu_h)\in {\mathit{\boldsymbol{V}}}^\ell_h\times W^\ell_h\times M^\ell_h$ and $i = 1, 2, \cdots, M$.

Proof. By the definition of $\Pi_V^\ell$ in Eq (4.14a), $\widehat \Pi^{\ell}$ in Eq (3.13b), and $\Pi_W^\ell$ in Eq (4.14b), we get

The initial identity has been established.

Going forward, we shall demonstrate the second identity. First, by the definition of $\Pi_W^\ell$ and $\Pi_V^\ell$ in Eq (4.14c), we have

and

Adding Eqs (4.17) and (4.18) yields

Therefore, combining with Eq (4.19), we have

This proves the second identity.

Subtracting the outcome of Lemma (4.3) from the FEMC-HDG system (3.12a)–(3.12b) yields the subsequent error equations.

Lemma 4.4. Let $\xi^{u^n}_{ih} = u^n_{ih} - \Pi^\ell_W u^n_i$, $\xi^{\mathit{\boldsymbol{p}}^n}_{ih} = \mathit{\boldsymbol{p}}^n_{ih} - \mathit{\boldsymbol{\Pi}}^\ell_V \mathit{\boldsymbol{p}}^n_i$, and $\xi^{\widehat{u}^n}_{ih} = \widehat{u}^n_{ih} - \widehat{\Pi}^\ell u^n_i$. Then we have that the following error equations:

and

hold for all $({\mathit{\boldsymbol{r}}}_h, v_h, \mu_h)\in {\mathit{\boldsymbol{V}}}^\ell_h\times W^\ell_h\times M^\ell_h$ and $i = 1, 2, \cdots, M$.

Here, we determine the estimates of $\mathcal{E}_{h}^{u^n}$, $\mathcal{E}_{M}^{u^n}$, $\mathcal{E}_{h}^{{\mathit{\boldsymbol{p}}}^n}$, and $\mathcal{E}_{M}^{{\mathit{\boldsymbol{p}}}^n}$ to establish an error assessment for the FEMC-HDG algorithm.

Theorem 4.2. Let $(u^n_i, \mathit{\boldsymbol{p}}^n_i)$ and $(u^n_{ih}, \mathit{\boldsymbol{p}}^n_{ih})$ denote the solutions to systems (1.2) and (3.12) at time $t_n$, respectively. Assume that $f_{i} \in \widetilde{L}^{2}(0, T; L^2(D))$, $g_{i} \in \widetilde{L}^{2}(0, T; L^2(\partial D_{1}))$, $u^0 \in \widetilde{L}^{2}(H^{\ell+2}(D))$, and that the condition $(\mathbf{H1})$ and $(\mathbf{H2})$ are met. Consequently, there exists a positive constant $C$ such that

provided the stability conditions $\kappa_{\min} - \kappa_{+} > 0$ and $\rho_{\min} - \rho_{+} > 0$ are satisfied.

Proof. First, by selecting $(\mathit{\boldsymbol{r}}_h, w_h, \mu_h) = (\xi^{\mathit{\boldsymbol{p}}^n}_{ih}, \xi^{u^n}_{ih}, \xi^{\widehat{u}^n}_{ih})$ in Eqs (4.20) and (4.21), we obtain

and

By summing Eqs (4.23) and (4.24), applying the polarization identity, and integrating over the probability space, we derive

Utilizing the Cauchy-Schwarz inequality, the conditions (3.6) and (3.4), and Young's inequality on the right-hand side of Eq (4.24), we obtain the following five inequalities for any $\epsilon_i > 0$ $(i = 1, 2, \cdots, 6)$:

and

By substituting Eqs (4.25) through (4.29) into Eq (4.24), and additionally applying the conditions (3.7), (3.8), (3.3), and (3.5), we obtain the following for any $0 < \sigma_1 < 1$:

According to the condition (4.1), we choose

Then, from Eq (4.30), we obtain

Additionally, the following estimates pertain to temporal errors:

and

Summing Eq (4.31) up from $n = 1$ to $N$ and using Eqs (4.32)–(4.35), we arrive at

Applying Gronwall's inequality to Eq (4.36), and utilizing Lemma 4.1 along with Lemma 4.2, we obtain

Decomposing $\mathcal{E}_{h}^{u^n}$ as

and $\mathcal{E}_{h}^{{\mathit{\boldsymbol{p}}}^n}$ as

and applying the triangle inequality, Jensen's inequality, and Lemma 4.2, we derive the expected result (4.22).

Remark 4.2. As far as we know, the previous other research has only provided a suboptimal $L^2$ convergence rate for ensemble solutions $u_{ih}$ for the heat equation with a Robin coefficient. In contrast, our result in Eq (4.22) achieves the optimal $L^\infty(0, T;L^2(D))$ convergence rate on a general polygonal domain $D$.

The statistical errors $\mathcal{E}_{M}^{u^n}$ and $\mathcal{E}_{M}^{{\mathit{\boldsymbol{p}}}^n}$ can be obtained using the standard estimation method (refer to ^[1]).

Theorem 4.3. Given that the condition $(\mathbf{H1})$, $(\mathbf{H2})$, and the stability condition (4.1) hold, along with $f_i \in \widetilde{L}^{2}(0, T; L^2(D))$, $g_i \in \widetilde{L}^{2}(0, T; L^{2}(\partial D_{1}))$, and $u^{0}_i \in \widetilde{L}^{2}(H^{\ell+2}(D))$, then there exists a positive constant $C$ such that

Proof. We first analyze $\mathbb{E}\left[\left\| \mathcal{E}_{M}^{u^n}\right\|^2_{\mathcal T_h}\right]$, for all $1 \leq n\leq N$. It is easy to see that

Given that $u_{h}^{n}(\cdot, \omega_{1}), \cdots, u_{h}^{n}(\cdot, \omega_{M})$ are i.i.d., the expectation of $\left(\mathbb{E}[u_{h}^{n}] - u_{ih}^{n}, \mathbb{E}[u_{h}^{n}] - u_{jh}^{n} \right)$ is zero when $i \neq j$. Thus, we obtain

Let ${\mathcal{M}} = u_{h}^{n}$ and $\bar{\mathcal{M}} = \mathbb{E}\left[ \mathcal{M} \right]$, from which we can infer

As a result, we arrive at

Regarding $\mathbb{E}\left[{\|\mathcal{E}^{{q\mathit{\boldsymbol{}}}^n}_M\|^2_{\mathcal T_h}}\right]$, the situation is similar. From Theorem 4.1, we derive the result (4.38).

By combining the space error, time error, and the MC sampling error, we can derive the total error of the FEMC-HDG algorithm.

Theorem 4.4. Let $f \in \widetilde{L}^{2}\left(0, T; L^2(D)\right)$, $g \in \widetilde{L}^{2}\left(0, T; L^{2}(\partial D_{1})\right)$, and $u^{0} \in \widetilde{L}^{2}\left(H^{\ell+2}(D)\right)$. Assume that the condition ($\mathbf{H1}$), ($\mathbf{H2}$), and the stability condition (4.1) are satisfied. Then there exists a constant $C > 0$ such that

Proof. Employing Young's inequality along with the triangle inequality on the first term of the left-hand side of Eq (4.40), we obtain

Applying Jensen's inequality to the first term on the right-hand side of the preceding inequality yields

Thus, by employing Theorem 4.1, Theorem 4.2, and Theorem 4.3, the desired outcome is obtained. Likewise, the term

on the left-hand side of Eq (4.40) can be deduced in a similar way.

Remark 4.3. The method outlined above attains first-order accuracy in time. Higher-order temporal accuracy can be achieved using numerical algorithms like the $\mathrm{BDF(k)}$ (backward differentiation formula of order $k$) scheme, with $k \geq 2$ ^[25,26,27], which can be adapted for use with ensemble algorithms for the uncertain heat equation. However, these may demand more stringent stability conditions than those outlined in Eq (4.1). A stable numerical scheme can be constructed using similar methods, and its convergence theory is also analogous.

4.3.   Postprocessing

As defined in ^[28], the element-by-element postprocessing is as follows: Find $u^{n*}_{ih} \in \mathbf P_{\ell+1}(K)$ such that for all

where

After such postprocessing, we can use $\frac{1}{M} \sum\limits^M_{i = 1} u_{ih}^{n*}$ in place of $\frac{1}{M} \sum\limits^M_{i = 1} u_{ih}^{n}$ in the FEMC-HDG algorithm. From the numerical experiments presented later, it can be observed that, after such postprocessing, the discrete solution achieves a super-convergent rate under certain conditions on the domain. For instance, a convex domain is adequate.

5.   Numerical test

In this section, we verify Theorem 4.4 through numerical simulations and highlight the advantages of the FEMC-HDG algorithm over the FNEMC-HDG approach. In particular, we examine the heat equation (1.1) with random coefficients, defined on the unit square $[0, 1]^2$, with boundary conditions $\partial D_{0} = \{x_2 = 0\} \cup \{x_2 = 1\}$ and $\partial D_{1} = \{x_1 = 0\} \cup \{x_1 = 1\}$. Our numerical simulations utilize the open-source software NGSolve ^[29], which can be accessed at https://ngsolve.org/.

We conduct the experiment using the exact solution

where $\omega$ is uniformly distributed over $[0, 1]$ and $t \in [0, 1]$. The coefficients are set as

5.1.   Tests of convergence order of time and space

The FEMC-HDG scheme is applied in this experiment to simulate the ensemble with $M = 30$. Define

where $(u^n_{ih, \mathcal{N}}, \mathit{\boldsymbol{p}}^n_{ih, \mathcal{N}})$ and $(u^n_{ih, \mathcal{E}}, \mathit{\boldsymbol{p}}^n_{ih, \mathcal{E}})$ denote the non-ensemble solution and ensemble solution, respectively. In the experiment, both time and space are partitioned uniformly. We use $\ell = 0$, where the time step is set to $\Delta t = h$. For $\ell = 1$, the time step is set to $\Delta t = h^3$. The spatial step size $h$ is varied from $\dfrac{1}{2^1}$ to $\dfrac{1}{2^5}$. The corresponding errors and convergence rates are shown in Table 1.

From Table 1 it can be seen that when $\Delta t = h$ and $\ell = 0$, the convergence rate of both the FEMC-HDG and FNEMC-HDG schemes is $\mathcal{O}(h + \Delta t) = \mathcal{O}(h)$. Additionally, when $\Delta t = h^3$ and $\ell = 1$, the convergence rate for both schemes is $\mathcal{O}(h^2 + \Delta t) = \mathcal{O}(h^2)$. Further observations reveal that, under the same spatial discretization parameter $h$, the errors for both the non-ensemble and ensemble schemes are of the same magnitude. Notably, after postprocessing, the convergence rate of $u^*_{ih}$ reaches $\mathcal{O}(h^3 + \Delta t) = \mathcal{O}(h^3)$ when $\Delta t = h^3$, which demonstrates spatial super-convergence.

We chose a sample size of $M = 100$, with a mesh size of $h = \dfrac{1}{2^6}$, $\ell = 1$, and a time step $\Delta t = \dfrac{1}{2^6}$. The mean solution at $t = 0.5$ is calculated. The results are displayed in Figure 1 (left). To assess the efficiency of the FEMC-HDG algorithm, we compare its results with those from simulations using the FNEMC-HDG algorithm, employing the same sample values. The difference in the mean solutions between FEMC-HDG and FNEMC-HDG is depicted in Figure 1 (right).

The difference between the FEMC-HDG and FNEMC-HDG algorithms is around $10^{-6}$, which indicates that the FEMC-HDG algorithm delivers a similar level of accuracy to that of the FNEMC-HDG algorithm. This suggests that the FEMC-HDG approach performs comparably in terms of precision.

5.2.   Test of reducing computation cost

For further evaluation of the FEMC-HDG algorithm, we adopt piecewise constant elements ($\ell = 0$). We conduct a comparative analysis with the FNEMC-HDG algorithm. For a fair comparison, results from both algorithms are obtained under the identical conditions that the same time step and spatial dimensions are adopted. As illustrated in Table 2, the evaluations are performed with a mesh size $\Delta t = \dfrac{1}{2^4}$, and a time size of $h = \dfrac{1}{2^4}$. Given the moderate size of the discrete system, LU decomposition was employed.

From Table 2, it is evident that the FEMC-HDG algorithm surpasses the FNEMC-HDG algorithm in terms of CPU time. The FEMC-HDG algorithm improves computational efficiency by approximately 70% compared to the FNEMC-HDG algorithm.

5.3.   Test of Monte Carlo convergence speed

Next, we investigate the convergence rate of MC for the FEMC-HDG algorithm. Using the FEMC-HDG algorithm, we calculate the mean of the solution with $M_0 = 12,000$ samples as the reference benchmark, where $\ell = 0$ and $h = \Delta t = \dfrac{1}{2^4}$.

By adjusting the values of $M$ in the FEMC-HDG simulations, we assess the approximation errors to the benchmark. Furthermore, this error analysis is conducted for $J = 10$ independent runs, and the mean of the resulting errors is calculated.

Let $U_{M, h}^{n, m}$ represent the FEMC-HDG solution at time $t_n$ for the $m$-th independent trial, given by

where $u_{ih}^{n, m}$ is the result of the $m$-th experiment using the FEMC-HDG scheme. We also define the error measure

We performed 10 independent runs and reported the experimental results of $\mathcal{E}^{\mathcal{MC}}$ for $M = 10, 30, 90,$ and $270$ in Figure 2. The convergence rate observed in the experiments with respect to $M$ is approximately $-0.5$, which is consistent with the results derived in Theorem 4.4.

6.   Conclusion

In this study, we introduced a fully discrete HDG ensemble Monte Carlo algorithm that effectively tackles Robin boundary conditions and stochastic diffusion in heat equations. In the preparatory phase, the coefficient matrix for the linear systems is computed and preserved once, albeit requiring substantial computational effort. This methodology facilitates the swift assembly of linear systems in the operational phase, independent of the ensemble size. This marks the inaugural application of such an algorithm to address heat problems characterized by random diffusion and Robin coefficients. The performance of this method has been thoroughly validated and assessed.

The model described by Eq (1.1) is also relevant in robust optimal boundary control challenges as state equations (see ^[30,31]). Efficient resolution of the model in Eq (1.1) is paramount for numerically addressing the associated optimal Robin boundary control issues. Both theoretical analyses and numerical tests confirm the efficacy of the FEMC-HDG algorithm in resolving (1.1). Additionally, the FEMC-HDG algorithm extends to nonlinear random heat equations, demanding further exploration into stability conditions. Therefore, it is pertinent to investigate additional applications of the FEMC-HDG algorithm in random contexts, which we highlight as a valuable direction for subsequent studies.

Author contributions

Jinjun Yong: Methodology, analysis and writing original draft; Changlun Ye: Software and analysis; Xianbing Luo: Review, editing, and supervision.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is partially supported by the National Natural Science Foundation of China (12461076) and funding from the Scientific Research Fund Project of Guizhou Education University (2024ZD007).

Conflict of interest

The authors declare there is no conflict of interest.