Analytical-Numerical Solution of a Parabolic Diffusion Equation Under Uncertainty Conditions Using DTM with Monte Carlo Simulations

A numerical method to solve a general random linear parabolic equation where the diffusion coefficient, source term, boundary and initial conditions include uncertainty, is developed. Diffusion equations arise in many fields of science and engineering, and, in many cases, there are uncertainties due to data that cannot be known, or due to errors in measurements and intrinsic variability. In order to model these uncertainties the corresponding parameters, diffusion coefficient, source term, boundary and initial conditions, are assumed to be random variables with certain probability distributions functions. The proposed method includes finite difference schemes on the space variable and the differential transformation method for the time. In addition, the Monte Carlo method is used to deal with the random variables. The accuracy of the hybrid method is investigated numerically using the closed form solution of the deterministic associated 1 Universidad de Los Andes, Mérida, Venezuela, gcarlos@ula.ve. 2 Universidad de Córdoba Montería, Colombia, aarenas@correo.unicordoba.edu.co. 3 Universidad EAFIT, Medellín, Colombia, mcogollo@eafit.edu.co . Universidad EAFIT 49| Analytical-Numerical Solution of a Parabolic Diffusion Equation Under Uncertainty Conditions Using DTM with Monte Carlo Simulations equation. Based on the numerical results, confidence intervals and expected mean values for the solution are obtained. Furthermore, with the proposed hybrid method numerical-analytical solutions are obtained.


Introduction
Differential equations, in general, describe the rate of change in the physical property of matter with respect to time and/or space.Real phenomena lead to complicated differential equations, which seldom have exact solutions.The difficulty of obtaining analytical solutions and the availability of fast computing power make numerical techniques attractive.However, numerical methods can have slow convergence and instability problems.
Mathematical models dealing with uncertainty in differential equations have been considered in recent decades in a wide variety of applied areas,

|50
Ingeniería y Ciencia such as physics, chemistry, biology, economics, sociology and medicine.Differential deterministic equations have been extensively studied, both from both analytical and numerical viewpoints.However, in many situations, equations with random inputs are better suited in describing the real behavior of quantities of interest than their counterpart deterministic equations.Randomness in the input may arise because of errors in the observed or measured data, variability in experiment and empirical conditions, uncertainties (variables that cannot be measured or missing data) or plainly because of lack of knowledge [1], [2].
Differential equations where some or all the coefficients are considered random variables or that incorporate stochastic effects (usually in the form of white noise) have been increasingly used in the last few decades to deal with errors and uncertainty and represent a growing field of great scientific interest [3], [4], [5], [6].Additionally, uncertainty can be considered on the initial conditions or source term.Applications of the aforementioned random differential equations are wave propagations in homogeneous media, systems and structures with parametric excitations, dynamics of imperfectly known systems in physics, epidemics, medicine and biology [7], [8], [4].
Analytical treatment of random differential equations has been done by [4].It is important to remark that in general the statistical moments such as mean and variance of the solution process cannot be determined easily, see [4,Ch. 6] for details.Several applications to real world problems that consider randomness or uncertainty have been developed [9], [10], [11] The present work combines finite difference schemes for the discretization space with the differential transformation method for the time discretization.In addition, randomness is introduced in the diffusion P DE which models several physical processes and to the best of our knowledge this whole process has not be done before.The differential transformation method has been applied in several works and recently has been extended successfully to random differential equations [12], [13].The randomness is incorporated since inaccuracies in the physical measurements can affect several inputs of the diffusion P DE such as diffusion coefficient, source term, boundary and initial conditions, and can thus introduce some degree of uncertainty.In real world applications, the probability density functions of those quantities can only be estimated from physical measurements.It is important to remark that even if these measurements are done with the utmost care, the measured values will differ somewhat and some statistical tests are necessary to find the probability density distributions of these measurements.
The most commonly used distribution is the Gaussian one [1].However, other distributions such as Uniform, Beta and Gamma are also used and may be more appropriate if enough data points are available to suggest those distributions.
The real world is much more complex than anything that can be created with arithmetic and logical operations [14, p.341 and p.342].Therefore it is necessary to use methods that include some real world complexity like randomness [15].The Monte Carlo method is one of the most versatile and widely used numerical methods to deal with randomness.The power of Monte Carlo simulation modeling allows us to include more complexity to the deterministic mathematical model by incorporating the impact of randomness on the dependent variables [16].Random effects and the variations produced by using different probability distributions can be studied using Monte Carlo simulations.Applications of the Monte Carlo method in different areas are given, for example, in [17], [18], [11].
The Monte Carlo method is used here with the aim of obtaining qualitative and quantitative behavior of the numerical solutions of the random diffusion P DE.The Monte Carlo simulation differs from traditional simulation in that the model input parameters are treated as random variables, rather than as fixed values.These parameters must be identified with their uncertainty ranges and shapes of their probability density functions prescribed.There are no restrictions on the shapes of the probability density functions, although most studies make use of the basic ones, such as normal, log-normal, uniform, or triangular.Maximum and minimum limits on each input parameter can be prescribed to prevent unrealistic selection of extreme values [19].
The main reason here to apply the Monte Carlo method is due to the simplicity and effectiveness with which it deal with random variables.The major challenge with the method is to efficiently carry out many realizations and then to summarize the results into a few useful values such as the expectation and higher moments of the solution stochastic processes [20].The Monte Carlo method for solving random differential equations can be

|52
Ingeniería y Ciencia described briefly as: • Generate sample values of the random input from their known or assumed probability density function.
• Solve the deterministic equation corresponding to each value.
• Calculate statistics, such as mean and variance, of the set of deterministic solutions.
The number of realizations or runs can be as large as desired.However, there is always a compromise that must be reached concerning the optimal number of runs, since each run takes a certain time on the computer.Different criteria can be used to choose the number of realizations and this is not a straightforward task.A useful criteria on to select the optimal number of realizations is to stop when the increasing them does not significantly change expected value and variance of the physical process solution.Other criteria include theoretical ones using Central Limit Theorem, t-Student distribution and other probabilistic tools.Finally, it is important to remark that the numerical simulation time can be easily reduced using parallel processes.
The paper is organized as follows.In Section 2 the diffusion P DE is presented.Section 3 is devoted to present the basic properties of the differential transformation method.Numerical results for the deterministic and random P DE, including studies of the accuracy, are presented in Section 4. Finally, in Section 5 discussion and conclusions are presented.

Diffusion differential equation
In this paper, we consider a general random linear parabolic equation where the diffusion coefficient, source term, boundary and initial conditions include uncertainty.The deterministic associated problem is the following parabolic equation: Physically speaking, this model describes the heat conduction procedure in a given inhomogeneous medium with some input heat source f (x, t).The coefficient p(x, t) represents heat conduction property namely the heat capacity.The type of boundary condition used occurs when the ends of the bar are at given temperatures.
In a limited set of cases, solutions can be obtained in closed form by separation of variables.Thus, the great importance of numerical methods to calculate approximate solutions.Numerical approximations should preserve the dynamical behavior of the exact solution of the model and care should be used to avoid spurious solutions and instabilities.It is well known that the relation ∆t ∆x 2 , where ∆t is the time step and ∆x the spatial grid size, is critical for the numerical stability of the explicit finite difference schemes for this class of P DE, which puts serious restrictions on the time step.Therefore it is necessary to develop schemes that are robust with a threshold of convergence greater than for the traditional schemes.
In this paper, we introduce a hybrid method that combines the finite difference numerical schemes and an analytic-numerical method, the differential transformation method, to solve the diffusion P DE (1i)-(1iv) with randomness.The differential transformation technique is used to transform the governing equations from the time domain into the spectrum domain, followed by use of the finite difference method to formulate discretized iteration equations appropriate for rapid computation.Unlike the traditional high-order Taylor series method, which requires a many symbolic computations, the present method involves simple iterative procedures in the spectrum domain.The approximate solution is then obtained from a the partial sum of the inverse process.

Basic definitions of differential transformation method
Pukhov [21], proposed the concept of differential transformation, where the image of a transformed function is computed by differential operations.It is different from traditional integral transforms such as Laplace and Fourier.This method becomes a numerical-analytic technique that formalizes the Taylor series in a totally different manner.The differential transformation is a computational method that can be used to solve linear and non-linear ordinary and partial differential equations with their corresponding initial and/or boundary conditions.A pioneer using this method to solve initial value problems is Zhou [22], who introduced it in a study of electrical circuits.Additionally, the differential transformation method (DT M ) has been applied to solve a variety of problems that are modeled with differential equations [23], [24], [25], [26], [27].Some authors mentioned that DTM can be seen as computer-specialized procedure to compute the Taylor series [28].
The method consists of, given a system of differential equations and related initial conditions, transforming these into a system of recurrence equations that finally leads to a system of algebraic equations whose solutions are the coefficients of a power series solution.The numerical solution of the system of differential equation in the time domain can be obtained in the form of a finite-term series in terms of a chosen system of basis functions.In this method, we take {t k } +∞ k=0 as a basis functions and therefore the solution is obtained in the form of a Taylor series.Other bases may be chosen, see [29].The advantage of this method is that it is not necessary to do linearization or perturbations.Furthermore, large computational work and round-off errors are avoided.It has been used to solve effectively, easily and accurately a large class of linear and nonlinear problems.However, to the best of our knowledge, the differential transformation has not been applied yet in seasonal epidemic models.For the sake of clarity, we present the main definitions of the DT M as follows: Definition 3.1.Let x(t) be analytic in the time domain D, then it has derivatives of all orders with respect to time t.Let where k belongs to a set of non-negative integers, denoted as the K domain.Therefore, (2) can be rewritten as ing.cienc., vol.11, no.22, pp.49-72, julio-diciembre.2015.
where X(k) is called the spectrum of x(t) at t = t i .
Definition 3.2.Suppose that x(t) is analytic in the time domain D, then it can be represented as Thus, the equation ( 4) represents the inverse transformation of X(k).
where k ∈ Z + ∪ {0}, then the function x(t) can be described as where M(k) = 0 and q(t) = 0. M(k) is the weighting factor and q(t) is regarded as a kernel corresponding to x(t).
Note, that if M(k) = 1 and q(t) = 1, then Eqs. ( 3), ( 4), ( 5) and ( 6) are equivalent.Definition 3.4.Let [0, H] be the interval of simulation with H the time horizon of interest.We take a partition of the interval [0, H] as {0 = t 0 , t 1 , ... , t n = H} such that t i < t i+1 and and its differential inverse transformation of X(k) is defined as follow

|56
Ingeniería y Ciencia For the interested readers, the operation properties of the differential transformation method can be found in [21], [22].
From the definitions above, we can see that the concept of differential transformation is based upon the Taylor series expansion.Note that, the original functions are denoted by lowercase letters and their transformed functions are indicated by uppercase letters.Thus, applying the DT M , a system of differential equations in the domain of interest can be transformed to an algebraic equation system in the K domain and, thus x(t) can be obtained by a finite-term Taylor series plus a remainder, i.e., where In many modeling situations, the computation interval [0, H] is not always small, and in order to accelerate the rate of convergence and to improve the accuracy of the calculations, it is necessary to divide the entire domain H into n subdomains.This process can be seen as a implementation of an analytic continuation process due to the fact that the range of convergence of the direct sum of the series is too limited [28].Moreover, in this case, the DTM is truly different from the raw Taylor series method which, stricto sensu, does not involve any consideration of analytic continuation [28].The main advantage of the domain splitting process is that only a few Taylor series terms are required to construct the solution in a small time interval H i , where H = n i=1 H i .It is important to remark that, H i can be chosen arbitrarily small if necessary.Thus, the system of differential equations can be solved in each subdomain [24].Considering the function x(t) in the first sub-domain (0 ≤ t ≤ t 1 , t 0 = 0), the one dimensional differential transformation is given by ing.cienc., vol.11, no.22, pp.49-72, julio-diciembre.2015.
Notice that since the series converges in the domain [ti, ti + 1], its sum provides the solution in this domain with a sensible accuracy.Moreover, if this is true for each value t i then one gets, step by step, an approximate solution in the whole domain [0, H], each initial value x(H i being (approximately) provided by the sum of the series at the second boundary of the preceding sub-domain.Therefore, the differential transformation and system dynamic equations can be solved for the first subdomain and X 0 can be solved entirely in the first subdomain.The end point of function x(t) in the first subdomain is x 1 , and the value of t is H 0 .Thus, x 1 (t) is obtained by the differential transformation method as Since that x 1 (H 0 ) represents the initial condition in the second subdomain, then X 1 (0) = x 1 (H 0 ).And so the function x(t) can be expressed in the second sub-domain as In general, in each i − 1 subdomain one gets, (13) Using the D spectra method described above, the function x(t) can be obtained throughout the entire domain.
4 Construction of the hybrid scheme using DT M and finite differences for reaction-diffusion P DE Our goal in this section is to construct a hybrid scheme for Eq.(1i) combining the DT M and finite differences as follows: we begin taking the differential transforms of both sides of the governing equations from the

|58
Ingeniería y Ciencia time domain into the spectrum domain, i.e., taking the differential transformation with respect to the time variable t only.Thus, from Eq. ( 1) one gets the following spectrum: where U(x, k), G 1 (k), G 2 (k) are the differential transform of u(x, t), g 1 (t), and g 2 (t), respectively.
The finite difference method is then applied to Eq. ( 14), which contains only derivatives with respect to the space coordinates x.The whole domain is divided into M equal subintervals of length ∆x of the interval 0 ≤ x ≤ 1.The x coordinates of the grid points are given by x j = j∆x, for j = 0, ..., M .Using the central difference formula on the first and second derivatives, and the convolution of transformation in Eq. ( 14), the corresponding difference equation is given by Thus, we obtain a numerical scheme to compute the coefficients of the power series to obtain the solution in the respective subinterval, which is given by: ing.cienc., vol.11, no.22, pp.49-72, julio-diciembre.2015.
U j (0) =g j , For j = 0, ..., M , where bold uppercase letters are for the transformed functions.Thus, from a process of inverse differential transformation, the solutions on each subdomain can be obtained using n + 1 terms for the power series, i.e., and then the solution is given by:

Numerical results
In this section, numerical comparisons between the hybrid method and the analytical solution are presented in order to show the accuracy of the hybrid method for solving the diffusion P DE.In addition, combining the hybrid method with the Monte Carlo simulations we investigate the effect of introducing randomness on the diffusion coefficient, source term, boundary and initial conditions.Numerical simulations are performed with different realizations, time step sizes and parameters of the probabilistic distributions.The expected solution is computed as the average of several random solutions.The numerical results are presented mostly at the value x = 0.5, since this is the middle point of the interval of interest and generally is the point where the numerical methods give less accurate results.Graphical results corroborate this last fact.The randomness is included assuming a uniform probabilistic distribution.However, the methodology it is easily extendable to other probabilistic distributions.It is important to mention that the methodology proposed here is three-folded.The first step is to

|60
Ingeniería y Ciencia obtain the numerical scheme to compute the coefficients of the power series.This step includes to obtain the spectra of the equation considered and its discretization.In the second step, the numerical solution for each realization is computed.Finally, based on the Monte Carlo method many realizations are obtained and some statistics regarding the ensemble are computed.Therefore, the computation time of the whole process depends on many variables such as, the spectra, time step (for DTM and the finite difference scheme) and number of realizations.However, it has been mentioned in other works that the DT M is faster than the multistage Adomian method, but the Runke-Kutta methods require less computation time in comparison with DTM and multistage Adomian [12].
2. A collection of subsets of Ω called σ-algebra F ( or Borel field) that satisfies: the empty set is an element of F, F is closed under complementation and under countable unions.
3. A probability function P with domain F, and γ is a random variable (r.v.) defined on the original sample space Ω with probability function.
A solution of Example (5.1) with randomness means that for each γ ∈ Ω, u(x, t, γ) is a solution of the deterministic problem obtained from Example (5.1) taking realizations of the involved r.v, and where derivatives and limits are regarded in the mean square sense, see [4] for details.Since the solutions are stochastic processes, we can rely on the Monte Carlo method to compute them.It is important to remark that here we consider randomness on the initial conditions, diffusion coefficient, source term and boundary conditions separately.
The Monte Carlo method (with Simple Random Sampling) for uncertainty analysis is quite simple.At first, the more important inputs have to be identified as shown in (18).Next, uncertainty ranges and shapes of their probability density functions need to be chosen.Here, we choose as a first approach to test the methodology with uniform probability density functions for each r.v.. Thus, Example (5.1) is solved for each value assigned to the random variables to obtain results using the probability density function prescribed.The process is repeated many times, with the values of the input random parameters chosen from the corresponding distributions, in order to obtain a large ensemble of solutions [19].In other words, we sample from the probability density function in order to assign that value to the random variable considered, and solve the Example (5.1) using the DTM for each number sampled.It is important to emphasize that the method can be used with the input parameters having an arbitrary shape of probability density function.
In order to compute the coefficients of the power series of the random solution in the algebraic system (15) the following spectra are introduced: where i denotes the i-th split domain.

Numerical comparisons in accuracy
In Table 1 we present the absolute errors at x = 0.5 comparing the hybrid method and the analytical solutions of the deterministic problem of Example (5.1) when a space step size ∆x = 0.05 is used.The accuracy of the DT M hybrid method is improved by using the h-refinement approach (reducing the time step size).
In Table 2 we present the absolute errors at x = 0.5 comparing the hybrid method and the analytical solutions of the deterministic problem of Example (5.1) using different space step sizes ∆x and a fixed time step size ∆t = 0.0001.The accuracy of the DT M hybrid method is improved in this case by using the h-refinement approach in the space variable x.However, as it is remarked in [31] high space segmentation may cause divergence phenomenon.On the other hand, the p-refinement approach (increasing the number of terms) does not increase the accuracy of the hybrid method for this specific reaction-diffusion P DE due to its structure.
On the left hand side of Figure 1 it can be seen that the hybrid method reproduces the correct behavior of the solution for the diffusion P DE.In addition, on the right hand side of Figure 1 it can be observed that the absolute value of the error of the hybrid method increases with time and the maximum values are around the middle of the interval [0, 1].Thus, the study of the error at x = 0.5 is justified.These results show the numerical consistency of the hybrid method.
Table 1: Numerical absolute errors using the hybrid method with 3-terms, different time step sizes and a space step size ∆x = 0.05 for the deterministic version of Example (5.1) at x = 0.5 Table 2: Numerical absolute errors using the hybrid method with 3-terms, different space step sizes and a fixed time step size h = 0.0001 for the deterministic version of Example (5.1) at x = 0.5 In addition, it can be seen in Figures 2 and 3 that the amplitude of the confidence intervals decreases and converges to the expected solution, as we increase the number of realizations.This fact is of paramount importance from a physical point of view since it means that no matter how large is the error of the initial conditions measure, the solution will approximate asymptotically to the deterministic solution.The right hand side of Figure 3 shows the expected solution and confidence intervals when it is assumed normal probabilistic distribution for the initial conditions.Notice that the numerical results are similar to the ones with uniform probabilistic distribution but with less dispersion as was expected from the Gaussian form of the normal distribution.

Randomness on the boundary condition
In this subsection some Monte Carlo simulations are performed in order to observe the impact of boundary conditions uncertainties on the solution of the random version of Example (5.1).

|66
Ingeniería y Ciencia In Figure 4 it can be seen that the expected solution has a similar behavior as the deterministic solution when a uniform probabilistic distribution is used for the boundary conditions.The right hand side of Figure 4 shows that the amplitude of the confidence interval does not increase with time since the variance of the solution is exactly the variance of γ.This is due to the fact that perturbing the boundary conditions leads to an exact solution u(x, t) = e x+t + γ.

Randomness on the source term
In Figure 5 it can be seen that the expected solution agrees with the deterministic solution when a uniform probabilistic distribution is used for the source term.Notice that the amplitude of the 95% confidence interval (i.e. he range in which 95% of the realizations are found) increases with time.Therefore, the uncertainty on the source term is propagated over the time, giving rise with large probability to unpredictable feasible solutions.This fact is important from the physical point of view since it means that careful attention must be paid to the measure or estimation of the source term in order to have realistic solutions.ing.cienc., vol.11, no.22, pp.49-72, julio-diciembre.2015.

Randomness on the diffusion coefficient
In this subsection some Monte Carlo simulations are performed in order to observe the impact of randomness of the diffusion coefficient on the solution of the random version of Example (5.1).In Figure 6 it can be seen that the expected solutions agree with the deterministic solution when an uniform probabilistic distribution is used for the diffusion coefficient.
Notice that the amplitude of the 95% confidence intervals increase with time.Therefore, the uncertainty on the diffusion coefficient is propagated over the time, giving rise with large probability to unpredictable feasible solutions.

Conclusions
In this paper we apply a hybrid numerical method to solve random general linear diffusion equations where the diffusion coefficient, source term, boundary and initial conditions include uncertainty.The hybrid method combines the finite difference schemes for the discretization in space and the differential transformation method for the time discretization.In order to obtain accurate numerical solutions it is necessary to consider three

|68
Ingeniería y Ciencia issues: the first is the time step size used in the differential transformation method.The second one is the step size in the space for the finite difference scheme and the last is the order of the differential method.The accuracy of the DT M hybrid method can be improved by using the h-refinement approach in time and space variables.In addition due to the structure of the considered P DE the p-refinement approach does not improve the accuracy of the solutions for more than 3-terms.Nevertheless, in general increasing the differential transform order gives more accurate solutions at the expense of more computation time.The diffusion P DE has been selected due to the fact that reactiondiffusion equations arise in many fields of science and engineering, and, in many cases, there are uncertainties due to data that cannot be known, or due to errors in measurements and intrinsic variability.In order to model these uncertainties some probability distributions functions are assumed for the diffusion coefficient, source term, boundary and initial conditions.
The effect of introducing randomness in the diffusion P DE is justified by the fact that the diffusion coefficient, source term, boundary and initial conditions have some degree of uncertainty.Therefore, the random diffusion P DE is investigated by means of the well known Monte Carlo method.Based on the numerical results, confidence intervals and expected mean values for the solution are obtained.These confidence intervals are proportional to the variance of the probabilistic distributions of the random variable assumed for the diffusion coefficient, source term, boundary and initial conditions.This means that the dynamics behavior of the diffusion physical process can be predicted with some probability despite the uncertainty of the diffusion coefficient, source term, boundary and initial conditions.Finally it is important to mention that Monte Carlo simulations give realistic values which are consistent with the results obtained for the deterministic diffusion P DE.Future works can be developed for more complex cases such a nonlinear system with two interacting scalar fields.

Figure 1 :|64 Ingeniería y Ciencia 5 . 2
Figure 1: Numerical solution for the deterministic version of Example (5.1) using 5 terms for the DT M with h = 0.001 and ∆x = 0.1 for the finite difference (left).Plotting of the error for different times and values of the space variable x. (right)

Figure 2 :
Figure 2: Numerical solution for the random version of Example (5.1) at x = 0.5 assuming that initial conditions are perturbed by a term that follows a uniform distribution [−1.0, 1.0] with 50 (left) and 200 (right) realizations.Solutions are computed using 5 terms in the DT M , h = 0.0001 and ∆x = 0.1 for the finite difference.

Figure 3 :
Figure 3: Numerical solution for the random version of Example (5.1) at x = 0.5 assuming that initial conditions are perturbed by a term that follows a uniform distribution [−1.0, 1.0](left) and a normal distribution [0, 1/1.96] (right) with 150 and 200 realizations respectively.Solutions are computed using 5 terms in the DT M , h = 0.0001 and ∆x = 0.1 for the finite difference.

Figure 4 :
Figure 4: Numerical solution for the random version of Example (5.1) at x = 0.5 assuming that boundary conditions are perturbed by a term that follows a uniform distribution [−0.1, 0.1] with 100(left) and 20(right) realizations.Solutions are computed using 5 terms in the DT M , h = 0.0001 and ∆x = 0.1 for the finite difference.

Figure 5 :
Figure 5: Numerical solution for the random version of Example (5.1) at x = 0.5 assuming that the source term has in addition an uncertainty term that follows an uniform distribution [−1.0, 1.0] and performing 100 realizations.It is computed using 5-term in the DT M , h = 0.0001 and ∆x = 0.1 for the finite difference.

Figure 6 :
Figure 6: Numerical solution for the random version of Example (5.1) at x = 0.5 assuming that the diffusion coefficient has in addition an uncertainty term that follows an uniform distribution [−1.0, 1.0] with 50(left) and uniform distribution [−0.1, 0.1] with 150(right) realizations.Solutions are computed using 5-term in the DT M , h = 0.0001 and ∆x = 0.1 for the finite difference.