1.
Introduction
Let (Ω,F,P) be a complete probability space, where Ω is a sample space, F⊂2Ω is the σ-algebra of events, and P:F→[0,1] is a probability measure. [0,T] is the temporal interval. D⊆R2 is an open, bounded, Lipschitz physical space. ∂D=∂D0⋃∂D1 is the union of disjoint portions for the boundary. The unit outward normal vector to ∂D is denoted by 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]×D×Ω⟶R conforming to equations almost surely ω∈Ω,
where Robin coefficients α, diffusion coefficients a, boundary condition u, source term f, and initial condition y0 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]. ˉα(t,x) and ˉa(t,x) stand for the ensemble mean for the Robin factors as well as the diffusion factors, respectively, and we denote the fluctuations, α′j:=ˉα−αj, a′j:=ˉa−aj. We adopt a uniform partition on the interval [0,T] as well as tn=nΔt,n=0,1,2,...,N. Let ynj,fnj,unj,anj, αnj, ˉan, ˉαn, a′nj and α′nj be the values of functions yj,fj,uj,aj, αj, ˉa, ˉα, a′j and α′j at t=tn, respectively. yj=y(t,x;ω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 yn+1j(j=1,2,…,J;n=0,1,…,N−1) such that for all v∈H1(ˉD),
After arrangement, that is
Denote a∗=max1⩽j⩽J‖a′jˉa‖∞ and α∗=max1⩽j⩽J‖α′jˉα‖∞. From [30], we base on the following numerical stability condition a∗<1 and α∗<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, dx,ds, and dt in some expressions will be omitted when there is no confusion. The boundaries ∂D0 and ∂D1 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 Δt, sample size J, or mesh size h.
2.1. Basic preliminary.
Let ‖⋅‖ and (⋅,⋅) be the L2(D) norm as well as inner product, respectively. Simultaneously, ‖⋅‖∂D and (⋅,⋅)∂D stand for the L2(∂D) norm and the inner product. The Sobolev space Ws,q(D) with the norm ‖v‖Ws,q(D), here s∈N (natural number set) and q≥1. We denote Hs(D)=Ws,2(D). Particularly, H1(D) is equipped with the norm ‖⋅‖1=‖⋅‖1,D, which is defined by
Let H−s(D) be the dual space of the bounded linear functions on Hs(D), with norm ‖f‖−s=sup0≠v∈Hs(D)(f,v)/‖v‖s. Denote
where |v|∞=esssupx∈D|v|.
(Ω,F,P) is a complete probability space. Z∈L1P(Ω) is a random variable. Denote
Let δ=(δ1,…,δd) is a d-tuple with the length |δ|=d∑i=1δi,δi∈N+. The stochastic Sobolev spaces ˜Ws,q(D)= LqP(Ω,Ws,q(D)) contains stochastic function v:Ω×D→R, which is measurable with respect to the product σ-algebra F⊗B(D). The norm of ˜Ws,q(D) is defined by
Let ˜Hs(D)=˜Ws,2(D)≃ L2P(Ω)⊗Hs(D).
We will use the tensor product Hilbert space
with its inner product
The induced norm is given by
Lemma 2.1. The norm ‖⋅‖1,∂D1 is defined by
This norm and the standard norm ‖⋅‖1 in H1(D) are equivalent. Particularly, there exist constants CP,C′P>0 such that ∀v∈H1(D)
Proof. cf. [14,15]. □
The weak solution of the problem (1.1) is defined as follows: a function y∈X is a weak solution of (1.1) if it satisfies the initial condition y(0,x,ω)=y0(x,ω)∈L2P(H1(D);Ω) and for T>0,
for all v∈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,x,ω)∈L2P(0,T;L2(D);Ω) is uniformly bounded and coercive, i.e., ∃amin, amax>0 such that Pro{ω∈Ω:a(t,x,ω)∈[amin,amax],∀(t,x)∈[0,T]×D}=1.
(A2) α=α(t,x,ω)∈L2P(0,T;L2(∂D1);Ω) is uniformly bounded and coercive too, i.e., ∃αmin, αmax>0 such that Pro{ω∈Ω:α(t,x,ω)∈[αmin,αmax],∀(t,x)∈[0,T]×∂D1}=1.
Denote the ensemble group of solution variables for the equation (1.1) by yj(t,x)=y(t,x;ωj), j=1,2,...,J, with corresponding boundary conditions. [0,T]= ⋃N−1n=0[tn,tn+1] is an isometric partition of the interval. These norms are utilized in the numerical error discussion: ∀−1≤s<∞,
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 , we let and , where , . Without causing confusion, we abuse the fluctuating component symbols and . Suppressing the spatial discretization, utilizing an implicit-explicit for time-direct to the model (1.1), we propose the VMEMC scheme: find such that
Using a standard finite element (FE) method to discretize the space, we obtain the following block linear system for each ensemble member :
where is the mass matrix, is the boundary mass matrix, and is the diffusion matrix.
The above linear system is equivalent to the following:
where is the resulting coefficient matrix (independent of ). The matrix is symmetric positive definite (SPD) since both , , and 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 is a quasi-uniform triangulation of the domain , such that . Let be the diameter of the element . Define . Let be the set of polynomials of degree . Denote the finite elements (FEs) space of the domain by
the FEs space of the trace space corresponding to by here the space is defined as the restriction of on , is the number of nodes on , is the number of nodes on . The spaces above satisfy the following approximation properties [2,5] : ,
or denotes the interpolation operator related to the spaces or .
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 such that
here the initial value which .
As a comparison, here let's list the full-discrete VAEMC scheme [30] : find such that,
here the initial value which .
By the way, the full-discrete variational Monte Carlo (VMC) scheme: find such that,
here the initial value which .
Remark 3.1. (1) The symbols and 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, , , 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., , . Let be the fully discrete approximate solution at time , , , , . 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:
Choose a group of random samples for the stochastic Robin coefficient, diffusion coefficient, source term, boundary, and initial conditions, , , , , and for per . Therefore, the solutions are i.i.d..
Denote , and . For and the j-th member, we find to satisfy the VMEMC scheme. One finds the FE solution after choosing an approximate FE space.
Take the VMEMC sample average as an approximation to the expectation . Given an object of interest , we are going to analyze these outputs for the variational ensemble systems, , , , 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 , , and hypothesises and are satisfied. Then the VMEMC method is unconditionally stable. Furthermore, the numerical solution of (3.1) holds
Proof. Taking 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 . Multiplying by and integrating over the probability space, we can yield
For Eq (4.5), dropping two positive terms and , useing estimate Eqs (4.6) and (4.7) with , we can obtain
Summing from to , 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 and , 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 and 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 . Noting that , we can divide the overall error as
Here the statistical error, , is estimated using the sample number . The second term, , 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 and .
Theorem 4.2. Assume , , and hypotheses and 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 is examined in the following. Recall, the true solution satisfies
The error at the time is denoted by , where is a Lagrange interpolation [2,5] of . Taking and subtracting Eq (3.1) from (4.10), respectively, we have the error equation:
here 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 in the following lemma.
Lemma 4.1. Suppose every 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 is defined by Eq (4.12). The first term of is estimated as
The second and third terms of are estimated as
The last two terms of are estimated as
Combining the above estimates Eqs (4.15)–(4.21), using and , 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 satisfies the assumptions of Lemma 4.1. Moreover, suppose is an approximation of with the accuracy of the interpolation. Then, the VMEMC solution of (3.1) satisfies
Proof. Noting , rearranging the equation (4.11) with , we have
Integrating over the probability space and applying the polarization identity, we arrive at
From Lemma 4.1 with , 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 and , using the estimate (4.25)-(4.31), noting that , , selecting for , in equation (4.24), then we can obtain,
Summing to , collecting the constants, and reorganizing the items, we obtain
Then,
Inequalities and 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 , , satisfies the assumptions of Lemma 4.1. Moreover, suppose is an approximation of 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 as the manufactured solution; here , is a random variable, , and . 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 . 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 , , . We consider . We use isometric time partition and linear FEs on uniform space partition. Moreover, we take from to with 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 from to with 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 and boundary function in this numerical test.
5.2. Example two
The transient heat model (1.1) with random factors is considered on . 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 are independent. is uniformly distributed on . are uniformly distributed in the interval . The Robin boundary condition, the source term, and initial conditions, are chosen to match the exact solution. Let
(1) Convergence orders for time and physic space. Analogy to Example One, we use linear FEs 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 () and Table 4 ().
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 , which is higher than the theoretical result (Theorem 4.4). This may be the cause of the higher regularity of source term and boundary function .
(2) Convergence order for MC sampling. We take the VMEMC solution choosing samples as our benchmark, vary the values of in the VMEMC simulations, and then evaluate the approximation errors based on the benchmark. We repeat such error analysis for independent replicas and compute the average of the output errors. Denote the VMEMC solution at time in the th independent replica by where is result of the VMEMC (3.1) in the th experiment. Denote
We run times and take the logarithm for the numerical results of at are depicted in the left subgraph of Figure 1. One group of values of 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 is close to , which coincides with our theoretical results in Theorem 4.4.
(3) Comparison. Using the mesh size and the sample size , we obtain the mean and the variance of the solutions at . 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 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 , 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 of the random coefficient as the implicit term, and take the fluctuation 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 is not very large, such as , 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 .
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.