Reliable Efficient Difference Methods for Random Heterogeneous Diffusion Reaction Models with a Finite Degree of Randomness

This paper deals with the search for reliable efficient finite difference methods for the numerical solution of random heterogeneous diffusion reaction models with a finite degree of randomness. Efficiency appeals to the computational challenge in the random framework that requires not only the approximating stochastic process solution but also its expectation and variance. After studying positivity and conditional random mean square stability, the computation of the expectation and variance of the approximating stochastic process is not performed directly but through using a set of sampling finite difference schemes coming out by taking realizations of the random scheme and using Monte Carlo technique. Thus, the storage accumulation of symbolic expressions collapsing the approach is avoided keeping reliability. Results are simulated and a procedure for the numerical computation is given.


Introduction
Dealing with deterministic partial differential equations (PDE), finite difference methods (FD) are probably the most used because they are easy to apply and fairly efficient, [1,2]. Trying to capture real world problems, the models introduced uncertainty in several ways, basically assuming that both data, parameters, initial and or boundary conditions are stochastic processes instead of deterministic functions, [3,4]. The uncertainty appears not only because of error measurements, but also considering heterogeneity of the media, material impurities, or even the lack of access to measurements [5,6]. Independently of the type of modelling the uncertainty, the consideration of partial differential equations models (PDEM) has particular challenges. In fact, it is necessary to compute not only the stochastic process solution or approximating stochastic process, but also their statistical moments, mainly the expectation and the variance. Integral transforms methods are efficient techniques to solve PDEM based on integration resources in fitting domains [7,8]. Another powerful technique suitable for models with complex geometries is the finite element method [9]. Iterative methods, for instance FD have particular troubles derived from the storage accumulation of intermediate levels when the computer manages symbolically the involved stochastic processes, [10][11][12]. This drawback of the iterative methods for solving PDEM occurs in both approaches, the one based on Itô calculus [13] the so-called stochastic differential approach (SDEA), as well as the one based on the mean square calculus [14] also called random differential equations approach (RDEA). To face this computational challenge we take a set of realizations of the model, then we solve each sampled problem using FD method. Finally, we use Monte Carlo technique, [15,16] to average the results of the deterministic sampled problems to compute the expectation and the variance of the approximating solution stochastic process. In our model the involved stochastic processes (s.p.'s) are defined in a complete probability space (Ω, F , P) and have n degrees of randomness [14] (p. 37), i.e., they only depend on a finite number n of random variables (r.v.'s): where A i , i = 1, . . . , n, are mutually independent r.v.'s; F is a differentiable real function of the variable s (being s the spatial variable x or the temporal one t).
In addition, under this hypothesis, the s.p. h(s) has sample differentiable trajectories (realizations), i.e., for a fixed event ω ∈ Ω, the real function h(s, ω) = F(s, A 1 (ω), A 2 (ω), . . . , A n (ω)) is a differentiable function of the real variable s. For the sake of clarity in the presentation and to save notational complexity, we will assume that involved s.p.'s in the coefficients and initial or boundary conditions, have one degree of randomness, i.e., they have the form with A a r.v. and F a differentiable real function of the variable s. Then the s.p. h(s) has sample differentiable trajectories, i.e., for a fixed event ω ∈ Ω, the real function h(s, ω) = F(s, A(ω)) is a differentiable function of the real variable s. This paper deals with random parabolic partial differential models of the form where the diffusion coefficient p(x), the reaction coefficient q(x), the boundary conditions g i (t), i = 1, 2, and the initial condition f (x) are s.p.'s with one degree of randomness in the sense as defined above. In addition we assume that p(x), q(x), f (x) and g i (t), i = 1, 2 are mean square continuous s.p.'s in variables x and t, respectively, p(x) is also a mean square differentiable s.p. and the sample realizations of the random inputs p(x), q(x), g i (t), i = 1, 2 and f (x) satisfy the following conditions denoting p (x) as the mean square derivative of p(x): |p This model is frequent in chemical engineering sciences and in heat and mass transfer theory for reaction-diffusion problems with parameters depending on the spatial variables as it occurs in heterogeneous and anisotropic solids, [17] (p. 455), [18] (p. 388), [19,20].
The paper is organized as follows. Section 2 deals with some preliminaries, definitions and notations about the mean square calculus as well as the construction of the random mean square finite difference scheme (RMSFDS) resulting from the discretization of model (3)- (6). Section 3 is addressed to the study of properties of the RMSFDS such as positivity, stability and consistency. Throughout a sample approach, sufficient conditions for stability and positivity of the random numerical solution s.p. in terms of the data and discretization step-sizes are found. Consistency of the RMSFDS with Equation (3) is also treated throughout a sample approach and the consistency of the corresponding realized deterministic scheme for each fixed event ω ∈ Ω. In Section 4 we construct an algorithm to perform the efficient computation of the expectation and the variance of the numerical solution s.p. using Monte Carlo method and the solution of the sampled scheme. Numerical simulations for a problem where the exact solution is known are performed showing the efficiency of the proposed numerical method as well as the computations of the expectation and the variance of the numerical approximated s.p. A conclusion Section 5 ends the paper.

Preliminaries and Construction of the Random Finite Difference Scheme
For the sake of clarity in the presentation, in this section we recall some definitions and concepts related to the L p -calculus, [14]. In a probability space (Ω, F , P ), we denote L p (Ω) the space of all real valued r.v.'s U : Ω → R of order p, endowed with the norm where E[·] denotes the expectation operator, f U the density function of the r.v. U and ω an event of sample space Ω. Given T ⊂ R, a family of t-indexed r.v.'s, say {V(t) : t ∈ T}, is called a stochastic process of order p (p-s.p.) if for each t ∈ T fixed, the r.v. V(t) ∈ L p (Ω). We say that a p-s.p.
Furthermore, if there exists a p-s.p. V (t), such that then we say that the s.p. {V(t) : t ∈ T} is p -th mean differentiable at t ∈ T and V (t) is the p-derivative of V(t). In the particular case that p = 2, L 2 (Ω), definitions above leads to the corresponding concept of mean square (m.s.) continuity and m.s. differentiability, respectively.
In this section, we construct an explicit random finite difference scheme for solving problem (3)- (6). Firstly, let us write Equation (3) into the following form where p(x) ∈ L p (Ω) is p-th mean continuous and differentiable, p (x) is the p-derivative of p(t) and q(x) ∈ L p (Ω) is p-th mean continuous. Let us consider the uniform partition of the spatial interval [0, 1], of the form x i = ih, 0 ≤ i ≤ M, with Mh = 1. For a fixed time horizon, T, we consider N + 1 time levels t n = nk, 0 ≤ n ≤ N with Nk = T. The numerical approximation of the solution s.p. of the random problem (3)-(6) is denoted by u n i , i.e., By using a forward first-order approximation of the time partial derivative and centred second-order approximations for the spatial partial derivatives in Equation (13) one gets the following random numerical scheme for the spatial internal mesh points where The resulting random discretized problem (3)-(6) can be rewritten in the following form Please note that all the inputs of the random problem (16) are s.p.'s depending on one finite degree of randomness and lie in L p (Ω).

Properties of the Random Numerical Scheme: Positivity, Stability and Consistency
We are going to prove the positivity of the random numerical solution u n i , 0 ≤ i ≤ M , 0 ≤ n ≤ N} of the random scheme (16) and its · p -stability in the sense of fixed station respect to the time. We extend this type of stability to the random field.

Definition 1.
A random numerical scheme is said to be · p -stable in the fixed station sense in the where C is independent of the step-sizes h, k and the time level n.
First, we are going to find sufficient conditions on the spatial step-size h and the temporal step-size k, so that the numerical solution {u n i (ω)} constructed by sampling random scheme (16) for a fixed ω ∈ Ω guarantee its positivity, i.e., u n i (ω) ≥ 0 for 0 ≤ i ≤ M and for each time level n, 0 ≤ n ≤ N. We denote then the sampling scheme (18) can be rewritten as follows To guarantee the positivity of the numerical approximation {u n i (ω)} it is sufficient to impose the positivity of coefficients defined in (19). Please note that the simultaneously positivity of coefficients Taking into account the bound condition (8)

it follows that coefficients
Please note that for the particular case where p i (ω) is constant the positivity of coefficients A i (h, k)(ω) and C i (h, k)(ω) defined in (19), is established for h > 0. To guarantee the positivity of coefficient B i (h, k)(ω) from (19) and bounds (7)-(9) note that Thus, the positivity of Then taking into account the sufficient conditions (21), (23) and (24) over the discretization step-sizes h and k, the positivity of all the coefficients (19) of sampling scheme (18) for a.e. ω ∈ Ω is guaranteed and consequently the positivity of the numerical solution Let us prove now that random numerical scheme (16) is · p -stable in the sense of Definition 1. In this study we need to distinguish two cases for the sampling parameter q i (ω) for a fixed ω ∈ Ω.

From (18) imposing conditions (21) and (24) one gets
where Using (11), the boundary conditions of (18) and (25) we have by recurrence We denote where Thus, from (27) and (28) we obtain the following upper bound for the numerical solution of sampling scheme (18) for each level n, 0 ≤ n ≤ N = Tk, and for each spatial point From (18) imposing conditions (21) and (23) and using (25) and (26) we obtain Then using the boundary conditions of (18) and applying recurrently the bound exhibits in (31) one gets Taking into account the following inequality (1 + k |q min |) s ≤ (1 + k |q min |) N < e T |q min | , 0 ≤ s ≤ N , and the notation introduced in (28), we obtain this upper bound for the numerical solution of sample scheme (18) 0 ≤ u n i (ω) ≤ e T |q min | G(T) , for a.e. ω ∈ Ω , for each level n, 0 ≤ n ≤ N = Tk, and for each spatial point Please note that both bounds (30) and (33) are independent of n, h and k.
Finally, under discretization step-size conditions (21), (23) and (24) and from the upper bounds (30) and (33) it follows that where G(T) is defined in (28) and (29) and Consequently, random numerical scheme (16) is · p -stable in the sense of Definition 1. Summarizing, the following result was established. To prove the consistency of the random finite difference scheme (16) with the random partial differential Equation (13) let us introduce the following definition inspired in the well-known concept of consistency for deterministic PDEs, see [2]. Definition 2. Let us consider a RMSFDS F(u n i ) = 0 for a random partial differential equation (RPDE) L(u) = 0 and let the local truncation error T n i (U(ω)) for a fixed event ω ∈ Ω be defined by T n i (U(ω)) = F(U n i (ω)) − L(U n i (ω)), where U n i (ω) denotes the theoretical solution of L(u)(ω) = 0 evaluated at (x i , t n ). We call T n i (U) by With previous notation, the RMSFDS F(u n i ) = 0 is said to be · p -consistent with the Next result shows the consistency in the p-norm of RFDS (16) with RPDE (13). (16) is · p -consistent with the RPDE (13).

Theorem 2. The RFDS
Proof. Please note that for each fixed ω ∈ Ω the local truncation error using (13) and (15) is given by Assuming that U(x, t)(ω) is four times continuously differentiable with respect to x and two times continuously differentiable with respect to t and using Taylor expansions of U(x, t)(ω) at (x i , t n ) one gets where t n < δ < t n+1 , x i−1 < η j < x i+1 , j = 1, 2. Let us denote As we are in the scenario of finite degree of randomness and the involved variables have a truncated range, there exist D j (i, n), j = 1, 2, 3, positive constants such that From Definition 2, condition (7) and (36)-(40) it follows that

Algorithm
From a computational point of view, as it was commented on in the Introduction Section, the handling of the random scheme (16) in a direct way makes unavailable the computation of approximations beyond a few first temporal levels. This is because, throughout the iterative temporal levels, n = 1, · · · , N, it is necessary to store the symbolic expressions of all the previous levels of the iteration process collecting big and complex random expressions with which the expectation and the standard deviation must be computed. Furthermore, although the random expressions can be stored it does not guarantee that the two first statistical moments could be computed in a numerical way. For this reason, we propose using the random numerical scheme (16) together with the Monte Carlo technique avoiding the described computational drawbacks. The procedure is as follows: to take a number K of realizations of the random data involved in the random PDE (3)-(6) according to their probability distributions; to compute the numerical solution, u n i (ω j ), j = 1, · · · , K, of the sampling deterministic difference schemes (18); to obtain the mean and the standard deviation of these K numerical solutions evaluated in the mesh points i = 1, · · · , M − 1, at the last time-level N, denoted respectively by Algorithm 1 summarizes the steps to compute the stable approximations of the expectation and the standard deviation of the solution s.p., u n i , generated by means of the sampling difference scheme (18) and the MC method.

Numerical Example
To illustrate the efficiency of our proposed method in this subsection we present a test example where the exact solution s.p. is available. We consider the problem (3)-(6) with the random coefficients and the following boundary and initial conditions that is, u(1, t) = e c t e 2 2 + a e t , t > 0, Hereinafter, we will assume that a and c are independent r.v.'s. Please note that p(x) in (44) is a s.p. with one degree of randomness verifying condition (2) and g i (t), i = 1, 2, in (45) are s.p.'s with two degree of randomness (because they involve both r.v.'s a and c) verifying condition (2). Furthermore all random input data p(x), q(x), g 1 (t), g 2 (t) and f (x) lie in L 2 (Ω) and they are m.s. continuous and p(x) is m.s. differentiable too. In addition, conditions (7)-(11) are satisfied with a 1 = 0.4 e −1 , a 2 = 0.6 e 0 , −0.55 ≤ q(x, ω) ≤ −0.45 , ω ∈ Ω , 0 ≤ f (x, ω) ≤ 3.69453 , ω ∈ Ω.
From [18] (Section 3.8.5.) the exact solution of problem (46)-(49) when both parameters a and c are deterministic, is given by In our context, both a and c are r.v.'s, and expression (50) must be interpreted as a s.p. Then, using the independence between r.v.'s a and c, the expectation and the standard deviation of s.p. (50) are given by being (53) Figure 1 shows the evolution of the expectation (51) and the standard deviation (52) and (53) of the exact solution s.p. (50) when both parameters a and c are considered as the r.v.'s described above.  Table 1 collects the RMSEs (Root Mean Square Errors) computed at the time instant T = Nk = 1 with the temporal step-size k = 0.0001 (N = 10,000) for M − 1 = 79 internal spatial points x i = ih, 1 ≤ i ≤ 79 with h = ∆x = 0.0125 in the domain [0.0125, 1], using the following expressions where E[u(x i , t N )] and Var[u(x i , t N )] are given by (51)-(53), respectively. The good behaviour of both approximations the expectation and the standard deviation as the K simulations increases was observed. That is, the accuracy of the approximations to both statistical moments increases when the number of MC simulations is growing. In this sense, Figure 2 reflects the improvement of the approximations considering the study of the relative errors computed by the expressions RelErr .
(57)  Table 2 shows the second complementary study, where we have fixed the number of MC simulations K, K = 1600, and we have refined the step-sizes h and k attending to the stability conditions (21) Tables 3 and 4. As a result, a good strategy to study the convergence of approximations consists of choosing step-sizes h and k verifying the stability conditions and take a big enough number of realizations K such that the error does not vary significantly when one increases the number of realizations. For problems with no available solution the error is changed by the deviation between two successive numerical solutions.    The use of MC method has allowed the obtainment of approximations to the expectation and the standard deviation of the solution s.p. u N i of the random difference scheme (16) at time horizon T = Nk for N not necessarily small. However, if we use the random numerical scheme (16) directly in this example with the step-sizes h = 0.05 (M = 20) and k = 0.002 verifying the stability conditions (21) and (23), troubles appear in the early time-level n = 3, that corresponds to time t n = 0.006. In this case the symbolic expressions for the random numerical solution u n i and (u n i ) 2 , for n = 3 are available and their correspond expectations too. However, the Mathematica© software can not compute the numerical expectation of u n i 2 for 2 ≤ i ≤ 18, n = 3, in consequence it is not possible to compute the approximation of the standard deviation for these internal points at t n = 0.006 and hence at no other later time.

Conclusions
The main target of this paper is to solve the challenge of storage accumulation and further computational breakdown dealing with FD methods for solving random PDEM. Our approach is based on a combination of Monte Carlo method and the solution of sampled deterministic methods using explicit FD schemes. Explicitness is necessary to compute the statistical moments of the approximate solution what disregards the implicit FD methods. We here use the explicit classic difference method, but the Crank-Nicolson semi-implicit approach could be tried, making an ad hoc analysis. Numerical analysis provides sufficient conditions for positivity, stability and consistency for the proposed RMSFDS. Numerical experiments illustrate the reliability of the approach.
Author Contributions: M.C.C., R.C. and L.J. contributed equally to this work. All authors have read and agreed to the published version of the manuscript.