An ensemble scheme for the numerical solution of a random transient heat equation with uncertain inputs

An ensemble-based time stepping scheme is applied to solving a transient heat equation with random Robin boundary and diffusion coefficients. By introducing two ensemble means of Robin boundary and diffusion coefficients, we propose a new ensemble Monte Carlo (EMC) scheme for the transient heat equation. The EMC scheme solves a single linear system including several right-side vectors at each time step. Stability analysis and error estimates are derived. Two numerical examples verify the theoretical results and the validity of the EMC method.


Introduction
We consider the numerical simulations to a transient heat problem with random diffusion coefficients and Robin coefficients: Seek satisfying almost surely, is a Lipschitz domain. The boundary is divided into two disjoint parts 0 and 1 , i.e., 0 1 . n is the unit outward normal vector to . stands for a complete probability space, where is the sample space, 2 is the -algebra of events, and 0 1 is a probability measure.
The problem (1) describes a random heat and mass transfer process, which comes from [3,4]. This process is affected by random initial temperature, random ambient temperature, random thermal conductivity (diffusion coefficient), and random convective heat transfer coefficient (Robin coefficient). This model (1) is a problem with random diffusion coefficient and Robin coefficient. The model (1) can also be found in other problems, such as [22,23].
For the problem with random diffusion coefficient, there are many numerical algorithms to solve them, such as polynomial chaos method, stochastic collocation method, and stochastic finite element method (see, e.g., [1,7,9,19,24,26,27]). Compared with the numerical methods to determine the problem, using these methods to solve the problem with random coefficients needs to rewrite complex programs. Monte Carlo (MC) method is another kind important method (see, e.g., [6,8,10,21]) for the problem with random coefficient. Once the random variables are sampled independently and identically distributed (i.i.d.), we can use the program of determine problem to solve them. After independent identically distributed sampling of random variables, problem (1) becomes a parameterized problem: Find such that . The calculation cost is very high since the number is relatively large. In order to improve computing efficiency, ensemble method is proposed and widely used (see, e.g., [11,16,20,21] and the reference therein ). The ensemble method change solving A z b 1 2 to solving Az b 1 2 . In [20], the authors studied the parabolic problem with random coefficients by using ensemble method, and obtained an error estimate. But the error estimate therein is not optimal with respect to (w.r.t.) space. In view of this problem, the authors of [17] combined the ensemble with HDG method to obtain an optimal error estimate in space. To improve efficiency, multi-level ensemble method and second-order scheme are adopted in [17,21]. For the heat conduction problem with random diffusion and Robin coefficient, we have not found relevant results with the EMC method to our best knowledge.
In this work, a variational ensemble Monte Carlo scheme is first applied to handle simultaneously random diffusion coefficient and Robin coefficient in the transient heat (1). To do this, we consider the numerical simulation of model (2) firstly because the numerical approximation of (1) can be obtained by averaging 1 2 ... . Denote the ensemble average for the diffusion coefficient and the Robin coefficient by Through a spatial discretization, we observe that the coefficient matrix of the resulting linear system will be independent of . The discrete systems share a common coefficient matrix, and the right-hand-side (RHS) vectors vary among the ensemble members. Then, if the scale of the problem is small, the solution of the group can be obtained by LU decomposition of the coefficient matrix only once (see, e.g., [20,21]). While the case is of a large-scale, block Krylov subspace iteration method will be used to compute efficiently (see, e.g., [5,14,25]).
The following describes the structure of this article. In Section 2, some notations and preliminaries are introduced. The full discretization ensemble scheme for (5) is given in Section 3 with its stability and convergence analysis. The random transient heat equation and its stability as well as error analysis are discussed in Section 4. In Section 5, two numerical tests are presented. Some conclusions are given in Section 6.

Basic preliminaries
In this section, we will give some notations. For simplicity, x s, and in some expression will be omitted when there is no confusion. The boundaries 0 and 1 concern to the experimentally accessible and inaccessible parts, respectively. Let and be the 2 norm as well as inner product, respectively. Simultaneously and stand for corresponding the 2 norm as well as inner product. The Sobolev space with the norm , here (positive integer set) and 1 . We denote 2 . Particularly, 1 is equipped the norm 1 1 , which is defined by Let is the dual space of bounded linear functions on , with norm sup 0 . The norm which is equivalent to the standard norm 1 , cf. [12,13]. In the later statements, it will not differentiate the norm ess sup x (see, e.g., [2,19] . We will use the Hilbert space The FEs space of the boundary by Throughout this work, is a positive constant, it has different values in different places, and does not rely on time step , sample size , and mesh size .

A variational ensemble scheme for deterministic transient heat equations
In this section, we first give some assumptions and a full discretization ensemble scheme for (5). Then the stability and error estimate results are presented.
Suppose the following two hypotheses H1 and H2 hold. H1 There exists positive constants , min and max , such that for 1 2 , H2 There exists positive constants , and , such that for 1 2 , and Evidently, hypothesis H1 means the problem is uniformly coercivity. Hypothesis H2 shows that the difference between x and x is uniformly bounded, and so does the Robin coefficient x . Under the assumption that the isometric time division on 0 , the full-discrete scheme for (5) is that: seek here , the initial value 0 , 0 0 .

Stability
The stability of the ensemble scheme (12) has the following result.

Remark 1
The stability condition (13) requests, for 1 , the difference between and its mean is smaller than the coercivity constant . Similar requirement holds for 1 . Without such case, one partitions the ensemble into smaller groups and applies the ensemble algorithm to each smaller groups, and in this process, one must maintain the stability condition in these smaller groups. here , , is the 2 projection of onto , namely, 0 for any . One thus applies the decomposition in (23) and gets 1 1 1 1 1 .
(24) Letting 1 , using the polarization identity, and coercivity (8) and (9) For the last three items on the RHS of (26), it is that Replacing by these inequalities in (26), using the uniform coercivity H1 and uniform bounded condition H2 , amounting from 0 to 1, and multiplying both sides of (26) by 2 , we obtain where we used the assumption that 0 0 , thus 0 0 0 similar formula holds on boundary 1 . Applying the regularity hypothesis as well as standard FE estimates of the 2 projection in norm 1 (see, e.g., Section 4.4 of [2]), that is, Since the triangle inequality, one gets inequality (21).
The ensemble solution of the transient heat equations about uncertain inputs is investigated in next section.

The ensemble scheme for the stochastic transient heat model
Applying the ensemble scheme to stochastic PDEs (1), we first take the MC approach to random sampling. After sampling with independent identically distributed (i.i.d.), we solve (2). Then, the solution for (1) is obtained from the expectation of the solution of (2). An ensemble-based Monte-Carlo algorithm (i.e., EMC algorithm) is proposed to quantifies uncertainty and improve its computational efficiency. This algorithm is comprised of three steps: Next, we give the stability analysis and error estimate.

Stability
Paralleling to handle the PDE problem (2), we select the same FE space and mentioned in Section 2. Let x . Given the jth sample and 0 1 1, the EMC finite element method is to seek an approximation solution 1 such that  (27) is a random variable. To analyze the corresponding stability and error estimate, we suppose the following two hypotheses H3 and H4 are satisfied.
H3 There exist three positive constants , min max , such that, for 0 , the probability The stability condition (32) confines the difference between diffusion coefficients and the mean. The Robin coefficient is similar. Parallel to the deterministic equation (see Remark 1), If the stability condition false, the ensemble can be divided into smaller groups to maintain the condition (32) that applies to each group, so the EMC scheme will apply to all the smaller groups.

Convergence analysis
The EMC approximate solution of the full discretization is defined as 1 1 . Thus, one will lead an estimate for in some averaged norms.
can be divide into two terms: here we apply . The first term is corresponding to the FE discretization error, , is controlled by the time step as well as mesh size.
For the second term that is the statistical error, , is controlled by the sample size . Next, we examine the bounds of and . And then we achieve an error estimate for the EMC scheme.

Theorem 4 Assume
is the solution to stochastic PDE (1) while with , suits to (27). The statistical error can be estimated through applying the standard error calculation of MC approach (see, e.g., [18]):  By Theorem 3, we obtain Combining the FE error and MC sampling error, we can get the error of EMC finite element method. (36) Proof For the first item on the LHS of (36), Young's inequality and triangle inequality imply Using Jensen inequality, one gets Thus, the result is obtained according to Theorems 4 and 5. Other items in the LHS of (36) be able to calculate by a similar way.

Numerical tests
Two numerical examples on the ensemble algorithm are listed. Example 1 is deterministic heat transfer model, which is intended to show Theorem 2. Example 2 is random transient heat equation, which is applied to verify Theorem 6 and reveals the effectiveness of the EMC Algorithm.

Example 1
The first experiment implements the ensemble scheme for the deterministic heat transfer model (2) and tests the numerical performance of the propose algorithm. 3. The exact solution is 1 cos 2 1 cos 2 2 exp 1 , where is a random perturbation in 0 1 , 0 1 and 1 2 0 1 2 . The diffusion coefficient and Robin coefficient are chosen as 2 1 sin sin 1 2 . The initial condition, Robin boundary function, and source term are selected to match the exact solution.
The ensemble scheme (12) is used to simulate the group in this experiment, which involves three members with 1 0.2, 2 0.7, and 3 0.8. Define We use linear FEs 1 and isometric time partition and uniform space partition. To check the convergence rate in time, we select the time step from 1 2 to 1 2 4 with 1 2 9 . The numerical results are listed in Table 1. To check the convergence order in space, we choose space step from 1 2 4 to 1 2 7 with 1 2 9 . The numerical errors are listed in Table 2.  Table 2, the rate of convergence w.r.t. space is about 2, which is higher than theoretical findings of Theorem 2. This may be the cause of higher regularity of source term and boundary function in this numerical test (In Theorem 2, 2 0 1 , 2 0 2 1 ). This numerical example shows that the ensemble method is effective for the parameterized problem (2) with diffusion coefficient and Robin coefficient.

Example 2
Here we consider transient heat equations (1) with random coefficients on unit square domain 0 1 2 . The exact solution, the diffusion coefficient, and Robin coefficient  Two realizations of 2 ... 5 are presented in Fig. 1. Figure 1 shows the uniform coercivity and uniform bound for Robin coefficient almost surely.
Analogy to deterministic case, we use linear FEs 1 and isometric time partition and uniform space partition. The corresponding numerical results to test the convergence rate of time and space are listed in Tables 3 and 4.  From Table 3, the convergence rate of time is agree with theoretical findings of Theorem 6. From Table 4, the convergence order of space is higher than our theoretical findings of Theorem 6. This also may be because the higher regularity of source term and part boundary function (In Theorem 6, 2 0 1 and 2 0 2 1 ). Next we check the convergence rate of sampling number for the EMC approximation. We choose the EMC solution using 0 5000 samples as our benchmark, vary the values of in the EMC simulations, and then evaluate the approximation errors based on the benchmark. Furthermore, we repeat this independent error calculation 10 times and compute the average of the output errors. Denote the EMC solution at time in the th independent replica by 1 1 , where is solution of the ensemble scheme (27) in the th experiment. Define We run 11 times, the numerical results of 2 at 10, 20, 40, 80, 160 are plotted in Fig. 2. It is seen that the rate of convergence with respect to is close to 0.5, which coincides with our theoretical results in Theorem 6.  Choosing the sample size 100, the mesh size 1 2 6 , 1 2 4 , we calculate the mean and variance of the solutions at 0.14. The results are presented in Fig. 3. To measure the effectiveness of the EMC scheme, we compare the result with that of individual finite element MC (IFE-MC) solutions using the same realizations. The difference between the mean and the IFE-MC solution is also shown in Fig. 3.
The difference is on the order of 10 3 , and the mean 1 . This implies that the EMC scheme can achieve basically accurate approximation as the IFE-MC implements.
To test the efficiency of our proposed ensemble scheme, we list the results of our EMC algorithm and MC method in Table 5 under the same time and space size. Here the discrete system's size is not too large, we apply MATLAB with its LU factorization.
From Table 5, we can see that our EMC algorithm requires less CPU time than MC method when 1. The EMC method improves the computational efficiency by about 50-60%, compared with the non-ensemble scheme. The advantage of EMC becomes more obvious when the ensemble size increases.

Conclusion
In this work, an ensemble MC method is firstly used to numerically solve a transient heat equation with random Robin and diffusion coefficient. One of the main contribution is that we propose a EMC finite element method to solve a stochastic transient heat equation with random Robin and diffusion coefficient. This method save the computation cost by introducing two ensemble means compared with MC method. The other contribution is that we obtain the stability analysis and error estimate of the EMC finite element method for the parameterized flow problem. The third contribution is that we obtain the stability analysis and error estimate of the EMC finite element method for the stochastic transient heat problem. The effectiveness of this method is tested by some numerical experiments. One disadvantage is that the experimental results on the convergence order of space are better than the theoretical results. This may be because our theoretical results are obtained under the condition that 2 0 1 , 2 0 2 1 In Theorem 2 and 2 0 1 , 2 0 2 1 . In Theorem 6, and the regularity of in the experiments are higher than that.
The EMC algorithm can be applied to other problem (such as stochastic convection diffusion problem) with random Robin and diffusion coefficient. It can also be used to solve stochastic Robin boundary control problem (e.g., [22,23]).
In the future work, combining the multi-level MC with ensemble method and higher order time scheme can be considered. Multi-level important sampling (see, e.g., [15]) with ensemble method is also another worth investigating.