Breast conserving surgery (BCS) with adjuvant radiation therapy showed improved prognosis compared with mastectomy for early staged triple negative breast cancer patients

1.   Introduction

Let $(\Omega, \mathcal{F}, P)$ be a complete probability space, where $\Omega$ is a sample space, $\mathcal{F} \subset 2^{\Omega}$ is the $\sigma$-algebra of events, and $P: \mathcal{F} \rightarrow [0, 1]$ is a probability measure. $[0, T]$ is the temporal interval. $D\subseteq \mathbb{R}^{2}$ is an open, bounded, Lipschitz physical space. $\partial D = \partial D_{0} \bigcup \partial D_{1}$ is the union of disjoint portions for the boundary. The unit outward normal vector to $\partial D$ is denoted by $\mathbf{n}$. We carry out the numerical simulations of a transient heat model with random Robin coefficients and diffusion coefficients. That's to find a random function $y:[0, T] \times D \times\Omega \longrightarrow \mathbb{R}$ conforming to equations almost surely $\omega \in \Omega$,

where Robin coefficients $\alpha$, diffusion coefficients $a$, boundary condition $u$, source term $f$, and initial condition $y^{0}$ are random fields with continuous and bounded covariance functions.

The model (1.1) describes the process of heat as well as mass transfer with random Robin factors and diffusion factors, which can be found in ^[3,4] and the references therein. These uncertain factors are involved in thermal environments. The temperature is quantified stochastically. In the review ^[3], this model can describe random high-cycle temperature fluctuations observed at the upper core structure of fast-breeder reactors or random variations in heat transfer coefficients around the stator vanes of gas turbines. The model (1.1) also appears in robust optimal boundary control problems as state equations, such as ^[25,26]. The effective implementation of these applications needs to solve (1.1) efficiently.

To solve PDEs with random coefficients, there are many numerical algorithms that include polynomial chaos method, stochastic collocation method, stochastic finite element method, MC method, and so on (see, e.g., ^{[1,9,11,21,22,27,29,31]}). Where the MC method (see, e.g., ^[8,10,12,24]) is easy to implement and the convergence property of the MC method is not dependent on the dimension of the parameters in a random model. When the MC method is adopted, we first need independent sampling and then deal with independent numerical simulation per sample. These simulations are affected by independent diffusion factors, Robin factors, forces, initial and boundary conditions. To improve calculating efficiency, the ensemble algorithms are widely applied (see, e.g., ^{[13,17,20,23,24]}). In ^[23,24], the authors studied the parabolic problem with random coefficients by using the ensemble method, and obtained an error estimate. But the error estimate therein is not optimal in physical space. The authors of ^[18] combined the ensemble MC with the HDG method to obtain an optimal error estimate in space. For the heat conduction model with random Robin factors, we have not found other relevant results about ensemble algorithms besides our recent work ^[30].

Look back at the numerical algorithm in the article ^[30]. $\bar{\alpha}(t, \mathbf{x})$ and $\bar{a}(t, \mathbf{x})$ stand for the ensemble mean for the Robin factors as well as the diffusion factors, respectively, and we denote the fluctuations, $\alpha_{j}^{\prime} : = \bar{\alpha}-\alpha_{j}$, $a_{j}^{\prime} : = \bar{a}-a_{j}$. We adopt a uniform partition on the interval $[0, T]$ as well as $t_{n} = n \Delta t, n = 0, 1, 2, ..., N$. Let $y_{j}^{n}, f_{j}^{n}, u_{j}^{n}, a_{j}^{n}$, $\alpha_{j}^{n}$, $\bar{a}^n$, $\bar{\alpha}^n$, $a_{j}^{\prime \; n}$ and $\alpha_{j}^{\prime \; n}$ be the values of functions $y_{j}, f_{j}, u_{j}, a_{j}$, $\alpha_{j}$, $\bar{a}$, $\bar{\alpha}$, $a_{j}^{\prime}$ and $\alpha_{j}^{\prime}$ at $t = t_{n}$, respectively. $y_j = y(t, \mathbf{x}; \omega_{j})$ and the others are similar. In ^[30], based on the implicit Euler method, we developed a variational average ensemble method for the random transient heat problem (1.1). Conveniently, we called the variational average ensemble Monte Carlo (VAEMC) method: find $y_{j}^{n+1}\; (j = 1, 2, \ldots, J; n = 0, 1, \ldots, N-1)$ such that for all $v \in H^{1}(\bar{D})$,

After arrangement, that is

Denote $a^* = \max \limits_{1 \leqslant j \leqslant J}\left\|\frac{a_j^{\prime}}{\bar{a}}\right\|_{\infty}$ and $\alpha^* = \max \limits_{1 \leqslant j \leqslant J}\left\|\frac{\alpha_j^{\prime}}{\bar{\alpha}}\right\|_{\infty}$. From ^[30], we base on the following numerical stability condition $a^* < 1$ and $\alpha^* < 1$ for the first-order ensemble time-stepping methods (1.2). In the ensemble algorithm in the literatures ^{[17,18,19,20,23,24,30]}, both stability and convergence were conditionally dependent on the ratio between the fluctuating and average values of the relevant factors.

Now, instead of employing ensemble means of Robin coefficients as well as diffusion coefficients in scheme (1.2), we propose a maximum ensemble of these coefficients at each time step. The process gets a shared coefficient matrix for the ensemble group to better improve computational efficiency. This is a new variational ensemble Monte Carlo numerical method, which we called the variational MAX ensemble Monte Carlo (VMEMC) method to conveniently distinguish it from other methods. In contrast with related methodologies, the novelty of this algorithm is that the numerical stability and convergence are unconditional.

The following is the outline of this work. Some mathematical preliminaries and notations are listed in Section 2. The full discretization VMEMC scheme for the random transient heat model (1.1) is given in Section 3. In Section 4, the unconditional numerical stability and error estimate are analyzed. Several numerical tests are shown in Section 5. Related discussions of this algorithm and corresponding summaries are presented in Section 6.

2.   Mathematical preliminaries

In this section, we will give some notations and mathematical preliminaries. For simplicity, $d\mathbf{x}, d\mathbf{s}$, and $d t$ in some expressions will be omitted when there is no confusion. The boundaries $\partial D_{0}$ and $\partial D_{1}$ concern the experimentally accessible and inaccessible parts, respectively. Throughout this work $C$ is a positive constant; it has different values in different places and does not rely on time step $\Delta t$, sample size $J$, or mesh size $h$.

2.1.   Basic preliminary.

Let $\|\cdot\|$ and $(\cdot, \cdot)$ be the $L^{2}(D)$ norm as well as inner product, respectively. Simultaneously, $\|\cdot\|_{\partial D}$ and $(\cdot, \cdot)_{\partial D}$ stand for the $L^{2}(\partial D)$ norm and the inner product. The Sobolev space $W^{s, q}(D)$ with the norm $\|v\|_{W^{s, q}(D)}$, here $s\in \mathbb{N}$ (natural number set) and $q \geq 1$. We denote $H^{s}(D) = W^{s, 2}(D)$. Particularly, $H^{1}(D)$ is equipped with the norm $\|\cdot\|_{1} = \|\cdot\|_{1, D}$, which is defined by

Let $H^{-s}(D)$ be the dual space of the bounded linear functions on $H^{s}(D)$, with norm $\|f\|_{-s} = \sup _{0 \neq v \in H^{s}(D)}(f, v) /\|v\|_{s}$. Denote

where $|v|_{\infty} = \operatorname{ess\; }{\sup }_{\mathbf{x} \in D}|v|$.

$(\Omega, \mathcal{F}, P)$ is a complete probability space. $Z\in L_{P}^{1}(\Omega)$ is a random variable. Denote

Let $\delta = \left(\delta_{1}, \ldots, \delta_{d}\right)$ is a $d$-tuple with the length $|\delta| = \sum \limits_{i = 1}^{d} \delta_{i}, \delta_{i}\in N^{+}$. The stochastic Sobolev spaces $\widetilde{W}^{s, q}(D) =$ $L_{P}^{q}\left(\Omega, W^{s, q}(D)\right)$ contains stochastic function $v: \Omega \times D \rightarrow R$, which is measurable with respect to the product $\sigma$-algebra $\mathcal{F} \otimes B(D)$. The norm of $\widetilde{W}^{s, q}(D)$ is defined by

Let $\widetilde{H}^{s}(D) = \widetilde{W}^{s, 2}(D) \simeq$ $L_{P}^{2}(\Omega) \otimes H^{s}(D)$.

We will use the tensor product Hilbert space

with its inner product

The induced norm is given by

Lemma 2.1. The norm $\|\cdot\|_{1, \partial D_{1}}$ is defined by

This norm and the standard norm $\|\cdot\|_{1}$ in $H^{1}(D)$ are equivalent. Particularly, there exist constants $C_{P}, C_{P}^{\prime} > 0$ such that $\forall v \in H^{1}(D)$

Proof. cf. ^[14,15]. □

The weak solution of the problem (1.1) is defined as follows: a function $y \in X$ is a weak solution of (1.1) if it satisfies the initial condition $y(0, \mathbf{x}, \omega) = y^{0}(\mathbf{x}, \omega) \in L_{P}^{2}\left(H^{1}(D); \Omega \right)$ and for $T > 0$,

for all $v \in X$.

We emphasize that the assumed regularity only requires the random fields to be square integrable. Assumptions on the finite-dimensional noise, that is, all the involved random input data depend on a finite number of real-valued random variables. Assumptions of input data for coefficients are considered:

(A1) $a = a(t, \mathbf{x}, \omega) \in L_P^2\left(0, T; L^2(D); \Omega\right)$ is uniformly bounded and coercive, i.e., $\exists \; a_{\min }$, $a_{\max} > 0$ such that $Pro\{\omega \in \Omega : a(t, \mathbf{x}, \omega) \in [a_{\min}, a_{\max}], \forall (t, \mathbf{x}) \in [0, T] \times D \} = 1$.

(A2) $\alpha = \alpha(t, \mathbf{x}, \omega) \in L_P^2\left(0, T; L^2\left(\partial D_1\right); \Omega\right)$ is uniformly bounded and coercive too, i.e., $\exists \; \alpha_{\min }$, $\alpha_{\max} > 0$ such that $Pro\{\omega \in \Omega : \alpha(t, \mathbf{x}, \omega) \in [\alpha_{\min}, \alpha_{\max}], \forall (t, \mathbf{x}) \in [0, T] \times \partial D_1 \} = 1$.

Denote the ensemble group of solution variables for the equation (1.1) by $y_j(t, \mathbf{x}) = y\left(t, \mathbf{x}; \omega_{j}\right)$, $j = 1, 2, ..., J$, with corresponding boundary conditions. $\left[0, T\right] =$ $\bigcup_{n = 0}^{N-1}\left[t_{n}, t_{n+1}\right]$ is an isometric partition of the interval. These norms are utilized in the numerical error discussion: $\forall -1 \leq s < \infty$,

We use the MC method for random sampling because it is non-intrusive, easy to implement, and its convergence is independent of the dimension of the uncertain model parameters. For $j\in\{1, 2, ..., J\}$, we let $a_{\max}^{n} = \max \limits_{1 \leq j\leq J} \sup\limits_{x\in\Omega} a_{j}\left(t_{n}, \mathbf{x}\right)$ and $\alpha _{\max }^{n} = \max \limits_{1\leq j\leq J}\sup \limits_{x \in \Omega}\alpha_{j}\left(t_{n}, \mathbf{x}\right)$, where $a_{j}\left(t, \mathbf{x}\right) = a\left(t, \mathbf{x}, \omega_{j}\right)$, $\alpha_{j}\left(t, \mathbf{x}\right) = \alpha\left(t, \mathbf{x}, \omega_{j}\right)$. Without causing confusion, we abuse the fluctuating component symbols $a_{j}^{\prime \; n}: = a_{\max }^{n}-a_{j}\left(t_{n}, \mathbf{x}\right)$ and $\alpha_{j}^{\prime \; n}: = \alpha_{\max }^{n}-a_{j}\left(t_{n}, \mathbf{x}\right)$. Suppressing the spatial discretization, utilizing an implicit-explicit for time-direct to the model (1.1), we propose the VMEMC scheme: find $y_{j}^{n+1} \in H^{1}(\bar{D})$ such that

Using a standard finite element (FE) method to discretize the space, we obtain the following block linear system for each ensemble member $j$ :

where $M$ is the mass matrix, $M_{\partial D_{1}}$ is the boundary mass matrix, and $D$ is the diffusion matrix.

The above linear system is equivalent to the following:

where $A$ is the resulting coefficient matrix (independent of $j$). The matrix $A$ is symmetric positive definite (SPD) since both $\frac{1}{\Delta t} M$, $a_{\max }^{n+1} D$, and $\alpha_{\max }^{n+1} M_{\partial D_{1}}$ are SPD. The system (2.4) can be solved with efficient block solvers ^[7,16,28]. Further, since only one coefficient matrix is required for computation each time step, the storage requirement is thereby reduced.

2.2.   Finite elements preliminary.

Suppose $\mathcal{T}_{h}$ is a quasi-uniform triangulation of the domain $\bar{D}$, such that $\bar{D} = \bigcup_{K \in \mathcal{T}_{h}} \bar{K}$. Let $h_{K}$ be the diameter of the element $K$. Define $h = \max _{K \in \mathcal{T}_{h}} h_{K}$. Let $\mathbf{P}_{k}$ be the set of polynomials of degree $k$. Denote the finite elements (FEs) space $V_{h}$ of the domain $\bar{D}$ by

the FEs space $S_{h}$ of the trace space corresponding to $V_{h}$ by $S_{h} = span\left\{\phi_{i}\right\}_{1}^{n_{\partial D_{1}}} \subset V_{h, \partial D_{1}},$ here the space $V_{h, \partial D_{1}}$ is defined as the restriction of $V_{h}$ on $\partial D_{1}$, $n_{\bar{D}}$ is the number of nodes on $\bar{D}$, $n_{\partial D_{1}}$ is the number of nodes on $\partial D_{1}$. The spaces above satisfy the following approximation properties ^[2,5] : $\forall 1 \leq l \leq k$,

$I_{h}: C(\bar{D}) \rightarrow V_{h}$ or $I_{h}: C(\partial D_{1}) \rightarrow S_{h}$ denotes the interpolation operator related to the spaces $V_{h}$ or $S_{h}$.

3.   variational MAX ensemble Monte Carlo (VMEMC) algorithm

We first present a full discretization VMEMC method with finite elements for the VMEMC scheme (2.2) of the random transient heat model (1.1). And then the VMEMC algorithm is proposed. The full-discrete VMEMC scheme: find $y_{j, h}^{n+1} \in V_{h}$ such that

here the initial value $y_{j, h}^{0} \in V_{h}$ which $\left(y_{j, h}^{0}, v_{h}\right) = \left(y_{j}^{0}, v_{h}\right), \forall v_{h} \in V_{h}$.

As a comparison, here let's list the full-discrete VAEMC scheme ^[30] : find $y_{j, h}^{n+1} \in V_{h}$ such that,

here the initial value $y_{j, h}^{0}\in V_{h}$ which $\left(y_{j, h}^{0}, v_{h}\right) = \left(y_{j}^{0}, v_{h}\right), \forall v_{h}\in V_{h}$.

By the way, the full-discrete variational Monte Carlo (VMC) scheme: find $y_{j, h}^{n+1}\in V_{h}$ such that,

here the initial value $y_{j, h}^{0} \in V_{h}$ which $\left(y_{j, h}^{0}, v_{h}\right) = \left(y_{j}^{0}, v_{h}\right), \forall v_{h} \in V_{h}$.

Remark 3.1. (1) The symbols $\alpha_{j}^{\prime n+1}$ and $a_{j}^{\prime n+1}$ both appear in VMEMC and VAEMC schemes. Although they all represent fluctuations in corresponding random coefficients, they all mean different things.

(2) Although the VMEMC and VAEMC schemes are formally similar, the novelty of the VMEMC scheme is its unconditional numerical stability and convergence. And the VMEMC scheme is more efficient in the numerical simulation.

(3) If the random thermal conductivity and convective coefficient are independent of time, such as, $\alpha: = \alpha(x, \omega)$, $a: = a(x, \omega)$, the accordingly analyzed analysis is analogous to that presented below. It (VMEMC) is unconditionally stable as well as the first-order accurate approach.

(4) Assume that the random thermal conductivity and convective coefficient are provided with uncertain temperature-dependent values as ^[6], e.g., $a: = a(y)$, $\alpha: = \alpha(y)$. Let $y_{h}^{n}$ be the fully discrete approximate solution at time $t_{n}$, $a_{h}^{n} = a\left(y_{h}^{n}\right)$, $a_{h}^{\prime \; n} = a_{\max }-a_{h}^{n}$, $\alpha_{h}^{n} = \alpha\left(y_{h}^{n}\right)$, $\alpha_{h}^{\prime \; n} = \alpha_{\max }-\alpha_{h}^{n}$. Then we can replace them accordingly in the VMEMC scheme. Similarly, unconditionally stable and first-order accurate methods under assumptions that two coefficients are continuously differentiable can also be obtained.

Applying the variational MAX ensemble scheme to the stochastic model (1.1), we first take the MC approach for random sampling. After sampling with independent identically distributed (i.i.d.), we solve the corresponding deterministic PDE for each sample. Furthermore, the solution for model (1.1) is made from the expectation of the solutions of the deterministic PDEs. A variational MAX ensemble-based Monte-Carlo algorithm (i.e., VMEMC algorithm) is proposed to quantifies uncertainty and improve its computational efficiency. The VMEMC algorithm includes the following three steps (S):

VMEMC algorithm:

$\mathbf{S1.}$ Choose a group of random samples for the stochastic Robin coefficient, diffusion coefficient, source term, boundary, and initial conditions, $\alpha_{j} \equiv \alpha(\cdot, \cdot, \omega_{j})$, $a_{j} \equiv a(\cdot, \cdot, \omega_{j})$, $f_{j} \equiv f(\cdot, \cdot, \omega_{j})$, $u_{j} \equiv u(\cdot, \cdot, \omega_{j})$, and $y_{j}^{0} \equiv y^{0}(\cdot, \omega_{j})$ for per $j\in \{1, 2, ..., J\}$. Therefore, the solutions $y(\cdot, \cdot, \omega_{j})$ are i.i.d..

$\mathbf{S2.}$ Denote $y_{j}^{n} = y(t_{n}, \mathbf{x}, \omega_{j})$, $a _{\max }^{n} = \max \limits_{1 \leq j \leq J} \sup \limits_{x \in \Omega} a_{j}(t_{n}, \mathbf{x})$ and $\alpha _{\max }^{n} = \max \limits_{1 \leq j \leq J} \sup \limits_{x \in \Omega} \alpha_{j}(t_{n}, \mathbf{x})$. For $n = 0, \ldots, N-1$ and the j-th member, we find $y_{j}^{n+1}$ to satisfy the VMEMC scheme. One finds the FE solution $y_{h}(\cdot, \cdot, \omega_{j})$ after choosing an approximate FE space.

$\mathbf{S3.}$ Take the VMEMC sample average $\frac{1}{J} \sum_{j = 1}^{J} y_{h}(\cdot, \cdot, \omega_{j})$ as an approximation to the expectation $\mathbb{E}[y]$. Given an object of interest $G(y)$, we are going to analyze these outputs for the variational ensemble systems, $G(y_{h}(\cdot, \cdot, \omega_{1}))$, $\ldots$, $G\left(y_{h}\left(\cdot, \cdot, \omega_{J}\right)\right)$, to obtain the stochastic information.

4.   Stability and error analysis

The unconditional stability and error estimates are presented in this section.

4.1.   Stability analysis

The unconditional stability for the VMEMC method is below.

Theorem 4.1. Assume $f_{j} \in$ $L^{2}\left(0, T; H^{-1}(D) \right)$, $u_{j}\in L^{2}\left(0, T; L^{2}(\partial D_{1}) \right)$, and hypothesises $(\mathbf{A1})$ and $(\mathbf{A2})$ are satisfied. Then the VMEMC method is unconditionally stable. Furthermore, the numerical solution of (3.1) holds

Proof. Taking $v_{h} = 2\Delta t y_{j, h}^{n+1}$ in (3.1), we obtain

Noting that

and

applying the polarization identity, integrating over the probability space, and rearranging terms in equation (4.2), we obtain

Combing Cauchy-Schwarz, Young's inequalities with the trace theorem and the Sobolev embedding theorem to the right-hand side (RHS) of (4.5), we have

Under Lemma 2.1, we obtain $\left\|\sqrt{a_{j}^{n+1}} \nabla y_{j, h}^{n+1}\right\|^{2}+ \left\|\sqrt{\alpha_{j}^{n+1}} y_{j, h}^{n+1}\right\|_{\partial D_{1}}^{2} \geq \frac{C_{P}^{2}}{2}\left\|\sqrt{a_{j}^{n+1}} y_{j, h}^{n+1}\right\|_{1}^{2}$. Multiplying by $\Delta t$ and integrating over the probability space, we can yield

For Eq (4.5), dropping two positive terms $\Delta t\mathbb{E} \left[\left\|\sqrt{a_{j}^{n+1}}\nabla y_{j, h}^{n}\right\|^{2}\right]$ and $\Delta t\mathbb{E} \left[\left\|\sqrt{\alpha_{j}^{n+1}} y_{j, h}^{n}\right\|_{\partial D_{1}}^{2}\right]$, useing estimate Eqs (4.6) and (4.7) with $\epsilon_{1} = \epsilon_{2} = \frac{C_{P}^{2}}{8}$, we can obtain

Summing from $n = 0$ to $n = N-1$, one obtains the stability result Eq (4.1). □

Remark 4.1. (1) The stability result above shows that we have controlled over the temperature approximation in both $L^{\infty}\left(0, T; L^{2}(\Omega)\right)$ and $L^{2}\left(0, T; H^{1}(\Omega)\right)$, unconditionally. Specifically, this numerical stability is independent of both the time step size and the grid size constraints, as well as the fluctuations in the random coefficients. This feature is different from the VAEMC approach ((see, e.g., ^[23,30]).

(2)Compared to the standard Backward Euler method, we see that the numerical dissipation is enhanced with the additional term

(3) If we consider the assumption of minimum regularity, that is,

hypotheses $(\mathbf{A1})$ and $(\mathbf{A2})$ are satisfied, then the VMEMC method also possesses unconditional stability.

4.2.   Error analysis

In this subsection, we first decompose the overall error of the VMEMC algorithm into the statistical error and the discretization error. Next, we sequentially estimate the statistical error (Theorem 4.2) and the discretization error (Theorem 4.3). Finally, we combine these two error estimates to obtain the overall error estimate of the VMEMC algorithm (Theorem 4.4).

The full-discrete VMEMC approximate solution of (1.1) is defined as $\Psi_{h}^{n} \equiv \frac{1}{J} \sum_{j = 1}^{J} y_{j, h}^{n}$. Noting that $\mathbb{E}\left[y^{n}\right] = \mathbb{E}\left[y_{j}^{n}\right]$, we can divide the overall error $\mathbb{E}\left[y^{n}\right]-\Psi_{h}^{n}$ as

Here the statistical error, $\mathcal{E}_{S}^{n} = \mathbb{E}\left[y_{j, h}^{n}\right]-\Psi_{h}^{n}$, is estimated using the sample number $J$. The second term, $\mathcal{E}_{h}^{n} = \mathbb{E}\left[y_{j}^{n}-y_{j, h}^{n}\right]$, that is corresponding to the total error of the time discretization and the FE discretization, is controlled by the time step as well as mesh size. Next, we examine the bounds of $\mathcal{E}_{S}^{n}$ and $\mathcal{E}_{h}^{n}$.

Theorem 4.2. Assume $f_{j} \in$ $L^{2}\left(0, T; H^{-1}(D) \right)$, $u_{j}\in L^{2}\left(0, T; L^{2}(\partial D_{1}) \right)$, and hypotheses $(\mathbf{A1})$ and $(\mathbf{A2})$ are held, for the numerical solution to the VMEMC (3.1), then

Proof. The process of this proof is similar to the Theorem 5 in ^[30]. □

The $\mathcal{E}_{h}^{n}$ is examined in the following. Recall, the true solution satisfies

The error at the time $t = t_{n}$ is denoted by $e_{j}^{n} = \left(y_{j}^{n}-I_{h} y_{j}^{n}\right)-\left(y_{j, h}^{n}-I_{h} y_{j}^{n}\right) = \phi_{j}^{n}-\psi_{j, h}^{n}$, where $I_{h} y_{j}^{n}$ is a Lagrange interpolation ^[2,5] of $y_{j}^{n}$. Taking $v = v_{h} \in V_{h}$ and subtracting Eq (3.1) from (4.10), respectively, we have the error equation:

here $\xi\left(y_{j}^{n+1}, v_{h}\right)$ is defined as

The following regularity assumptions are needed:

As a preparatory step for estimating the discretization error, we first estimate each term of the $\xi\left(y_{j}^{n+1}, v_{h}\right)$ in the following lemma.

Lemma 4.1. Suppose every $y_{j}$ fulfills the equation (1.1) and the regularity assumptions (4.13), then,

Proof. Applying Lemma 2.1, the Cauchy-Schwarz and Young's inequalities, the Taylor's Theorem with integral remainder, we obtain the following results for $\xi\left(y_{j}^{n+1}, v_{h}\right)$ is defined by Eq (4.12). The first term of $\xi\left(y_{j}^{n+1}, v_{h}\right)$ is estimated as

The second and third terms of $\xi\left(y_{j}^{n+1}, v_{h}\right)$ are estimated as

The last two terms of $\xi\left(y_{j}^{n+1}, v_{h}\right)$ are estimated as

Combining the above estimates Eqs (4.15)–(4.21), using $\left\|\nabla v_{h}\right\| \leq\left\|v_{h}\right\|_{1}$ and $\left\| v_{h}\right\|_{\partial D_{1}} \leq\left\|v_{h}\right\|_{1}$, and integrating over the probability space, we can yield Eq (4.13). □

With the results of Lemma 4.1, the major convergence result of this work can be proved next.

Theorem 4.3. Suppose $y_{j}$ satisfies the assumptions of Lemma 4.1. Moreover, suppose $y_{j, h}^{0} \in V_{h}$ is an approximation of $y_{j}^{0}$ with the accuracy of the interpolation. Then, the VMEMC solution of (3.1) satisfies

Proof. Noting $e_{j}^{n} = \phi_{j}^{n}-\psi_{j, h}^{n}$, rearranging the equation (4.11) with $v_{h} = 2 \Delta t \psi_{j, h}^{n+1}$, we have

Integrating over the probability space and applying the polarization identity, we arrive at

From Lemma 4.1 with $v_{h} = \psi_{j, h}^{n+1}$, the first term on the RHS of Eq (4.24) can be expressed as

Applying Lemma 2.1, the Taylor's theorem with integral remainder, the Young's inequality, and the Cauchy-Schwarz inequality on the second term of RHS of (4.24), we have

Using the Young's inequality and the Cauchy-Schwarz inequality to the third and fourth terms, we obtain

Using the Young's inequality and the Cauchy-Schwarz inequality as well as the trace theorem to the two boundary terms, we have

Under Lemma 2.1, proceeding in similar fashion with equation (4.8), we obtain

Dropping the positive terms $\Delta t\mathbb{E} \left[\left\|\sqrt{a_{j}^{n+1}}\nabla \psi_{j, h}^{n}\right\|^{2}\right]$ and $\Delta t\mathbb{E} \left[\left\|\sqrt{\alpha_{j}^{n+1}} \psi_{j, h}^{n}\right\|_{\partial D_{1}}^{2}\right]$, using the estimate (4.25)-(4.31), noting that $\|\nabla v\| \leq\|v\|_{1}$, $\| v\| \leq\|v\|_{1}$, selecting $\epsilon_{i} = C_{P}^{2} / 40$ for $i = 3, 4, ..., 12$, in equation (4.24), then we can obtain,

Summing $n = 0$ to $n = N-1$, collecting the constants, and reorganizing the items, we obtain

Then,

Inequalities $\left\|e_{j }^{n}\right\|^{2} \leq 2\left\|\phi_{j }^{n}\right\|^{2}+2\left\|\psi_{j, h}^{n}\right\|^{2}$ and $\left\|\phi_{j}^{n+1}-\phi_{j}^{n}\right\|^{2}\leq C \Delta t \int_{t_n}^{t_{n+1}} \left\|\phi_{j, t}(s)\right\|^{2} d s$ imply the result (4.22). □

Combining the MC sampling error and the discretization error, one can acquire the total error of the VMEMC FE approximation.

Theorem 4.4. Assume that $f_{j} \in$ $L^{2}\left(0, T; H^{-1}(D) \right)$, $u_{j}\in$ $L^{2}\left(0, T; L^{2}(\partial D_{1})\right)$, $y_{j}$ satisfies the assumptions of Lemma 4.1. Moreover, suppose $y_{j, h}^{0} \in V_{h}$ is an approximation of $y_{j}^{0}$ with the accuracy of the interpolation. Then,

Proof. In view of the first term on the left-hand side (LHS) of (4.35), the triangle inequality and Young inequality suggest

Using Jensen inequality to the first term on the RHS of the above inequality, one gets

Thus, the result can be obtained by Theorem 4.2 and Theorem 4.3. The remaining terms of the LHS of (4.35) can be handled similarly. □

Remark 4.2. Collecting these constants, the RHS of (4.35) can be expressed as

5.   Numerical tests

For comparison purposes, here we use the numerical examples in our recent work ^[30] to test the various performances of the VMEMC algorithm. Example One is a deterministic heat transfer model with a known exact solution, which is intended to show the convergence rate of time and space, respectively. A random transient heat equation is in Example Two, which is used to test Theorem 4.4 and reveal the effectiveness of the VMEMC algorithm.

5.1.   Example one

We take $y_j = (1+\omega_j)\exp(t+1)\cos(2\pi x_1)\cos(2\pi x_2)$ as the manufactured solution; here $j \in \{1, 2, 3\}$, $\omega_j\in [0, 1]$ is a random variable, $(x_1, x_2)\in[0, 1]^2$, and $t \in [0, 1]$. The upper and lower boundary of the domain are zero Neumann boundary; the left and right boundary are Robin boundaries. The Robin coefficient and thermal conductivity diffusion factor are specified by $\alpha_j = a_j = 2+(1+\omega_j)\sin(x_1x_2)\sin(t)$. The source term, initial condition, and Robin boundary function should be adjusted appropriately.

In this experiment, the ensemble schemes (3.1) and (3.2) are applied to simulate the ensemble set that includes three samples where $\omega_1 = 0.2$, $\omega_2 = 0.7$, $\omega_3 = 0.8$. We consider $\mathcal{E}^{E, j}_{L^2} = \boldsymbol{\||} y_j^n-y^n_{j, h} \boldsymbol{\||}_{2, 0}, j = 1, 2, 3$. We use isometric time partition and linear FEs $(k = 1)$ on uniform space partition. Moreover, we take $\Delta t$ from $1/2$ to $1/2^{4}$ with $h = 1/2^{9}$ to verify the convergence order in the temporal direction. The numerical results on the VMEMC and the VAEMC schemes are shown in Table 1. Meanwhile, we select $h$ from $1/2^{4}$ to $1/2^{7}$ with $\Delta t = 1/2^{9}$ to consider the convergence order in space. The resulted numerical errors on the VMEMC and the VAEMC schemes are listed in Table 2.

Both the VMEMC and the VAEMC algorithms, the convergence orders on temporal direction are about one from Table 1, and coincide with theoretical findings. From Table 2, the rates of convergence on the VMEMC scheme in spatial direction are close to 2, which is higher than the theoretical result (Theorem 4.4). This may be the cause of the higher regularity of source term $f$ and boundary function $u$ in this numerical test.

5.2.   Example two

The transient heat model (1.1) with random factors is considered on $[0, 1]^2$. The upper and lower boundary of the domain are zero Neumann boundary; the left and right boundary are Robin boundaries. Inspired by numerical examples in literature ^[25,30], we choose the exact solution, the diffusion factor and Robin factor are as follow:

The random variables $\omega_1, \cdots, \omega_5$ are independent. $\omega_1$ is uniformly distributed on $[0, 1]$. $\omega_2, \cdots, \omega_5$ are uniformly distributed in the interval $[-0.5, 0.5]$. The Robin boundary condition, the source term, and initial conditions, are chosen to match the exact solution. Let $\mathcal{E}_{L^2} = \sqrt{\frac{1}{J}\sum_{j = 1}^{J}\boldsymbol{\||} y_j^n-y^n_{j, h} \boldsymbol{\||}^{2}_{2, 0}}, j = 1, 2, 3.$

(1) Convergence orders for time and physic space. Analogy to Example One, we use linear FEs $(k = 1)$ on uniform partition physical space and isometric partition in time direction. The corresponding average numerical errors on VAEMC and VMEMC schemes are shown in Table 3 ($h = \frac{1}{2^8}, J = 50$) and Table 4 ($\Delta t = \frac{1}{2^9}, J = 50$).

From Table 3, the convergence order on temporal direction is according to Theorem 4.4. Looking at Table 4, the orders of convergence on the VMEMC scheme on spatial direction are close to $2$, which is higher than the theoretical result (Theorem 4.4). This may be the cause of the higher regularity of source term $f$ and boundary function $u$.

(2) Convergence order for MC sampling. We take the VMEMC solution choosing $J_{0} = 5000$ samples as our benchmark, vary the values of $J$ in the VMEMC simulations, and then evaluate the approximation errors based on the benchmark. We repeat such error analysis for $R = 12$ independent replicas and compute the average of the output errors. Denote the VMEMC solution at time $t_n$ in the $r-$th independent replica by $\Psi_{J, h}^{n, r} \equiv \frac{1}{J} \sum_{j = 1}^{J} y_{j, h}^{n, r},$ where $y_{j, h}^{n, r}$ is result of the VMEMC (3.1) in the $r-$th experiment. Denote

We run $R = 12$ times and take the logarithm for the numerical results of $\mathcal{E}^{*}_{L^2}$ at $J = 10, 20, 40, 80,160$ are depicted in the left subgraph of Figure 1. One group of values of $\mathcal{E}^{*}_{L^2}$ can be seen in Table 5 and the corresponding linear regression model is shown in the right subgraph of Figure 1.

We apply linear regression analysis on the numerical errors, which show

These results shows that the rate of convergence with respect to $J$ is close to $-0.5$, which coincides with our theoretical results in Theorem 4.4.

(3) Comparison. Using the mesh size $h = \Delta t = 1/2^{6}$ and the sample size $J = 100$, we obtain the mean and the variance of the solutions at $t = 0.2$. The results are plotted in Figure 2 and Figure 3 for different algorithms.

From Figure 2 and Figure 3, it is seen that both VMEMC and VAEMC approximations are excellent.

To check out the performance of the VMEMC algorithm, we compare the result with that of individual finite element MC (IFE-MC) simulations applying the same groups of sample values. The difference from the mean to the IFE-MC solution is plotted in Figure 4.

The accurate approximation to the IFE-MC implemented for the VMEMC scheme is very nice. From Figure 4, the difference on the order of $10^{-3}$ signifies that the VMEMC approximation can also achieve basically accurate approximation as the IFE-MC implements.

(4) Efficiency. To test the efficiency, setting the same time as well as space size, we record the CPU running times of the VMEMC, VAEMC, and MC algorithms in Table 6. Here we apply LU factorization because the discrete system's size is small.

From Table 6, when $J > 1$, we can show that both the VMEMC and the VAEMC algorithms take fewer CPU times than the MC algorithm. In contrast with the non-ensemble scheme, these ensemble methods improve the calculational efficiency by about 70%–80%. As the size of the ensemble increases, the efficiency of these ensemble algorithms improves, particularly for the VMEMC algorithm.

6.   Discussion and conclusion

The core idea of the ensemble algorithm is to develop a semi-implicit and semi-explicit temporal discretization scheme for PDEs. It will select a representative $\hat{a}$ of the random coefficient as the implicit term, and take the fluctuation $\| a_{j}-\hat{a}\|$ between random coefficients and this representation as the explicit term. This results in a shared coefficient matrix for each ensemble member at every time step. Consequently, storage requirements are reduced, and solution speed is increased. In general, as long as the relative fluctuation $\frac{ \| a_{j}-\hat{a}\|}{ \hat{a}}$ is not very large, such as $\frac{ \| a_{j}-\hat{a}\|}{ \hat{a}} < 1$, the ensemble numerical scheme can be stable. A large amount of literature ^{[17,18,19,20,23,24,30]} discussed the numerical stability conditions and error analysis of the ensemble algorithm based on taking the mean as representation $\hat{a}$.

Based on taking the maximum value of the random coefficient as representation, an unconditionally stable ensemble scheme, VMEMC, is established in this work. First, we apply the VMEMC algorithm to deal with the transient heat equation with the stochastic diffusion factors and the Robin factors. Without any stability conditions, we prove that the algorithm is numerically stable. Additionally, we also provide the overall error estimate for the algorithm. Moreover, numerical experiments are conducted to demonstrate the application of the ensemble schemes and to validate the established properties of the algorithm.

Important next steps include accommodating phase changes in the solid material (e.g., transitions to a liquid phase) and incorporating additional physical processes into the boundary conditions (e.g., surface-to-ambient radiation). The stochastic Robin boundary control problem is worth investigating. The VMEMC method can also be developed for nonlinear stochastic parabolic equations.

Author contributions

Tingfu Yao: Methodology, writing original draft, software; Changlun Ye: Software and supervision; Xianbing Luo: Review, editing, and supervision; Shuwen Xiang: Supervision and validation.

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 supported by the National Natural Science Foundation of China (Granted No. 11961008, 12361033) and The Natural Science Research Foundation of Education Department of Guizhou Province (Granted NO. QJJ[2024]190).

Conflict of interest

The authors declare there is no conflict of interest.