# Solucion numerico-analitica de una ecuacion de difusion bajo condiciones de incertidumbre utilizando DTM y Monte Carlo.

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

1 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, 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 PDE 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 PDE 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 PDE. 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 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 PDE is presented. Section 3 is devoted to present the basic properties of the differential transformation method. Numerical results for the deterministic and random PDE, including studies of the accuracy, are presented in Section 4. Finally, in Section 5 discussion and conclusions are presented.

2 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:

[u.sub.t](x, t) = [[p(x, t)[u.sub.x](x, t)].sub.x] + f (x, t), 0 < x < 1, t > 0, (1i)

u(0, t) = [g.sub.1](t), t > 0, (1ii)

u(1,t) = [g.sub.2](t), t> 0, (1iii)

u(x, 0) = g(x), 0 [less than or equal to] x [greater than or equal to] 1. (1iv)

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 [DELTA]t/[DELTA][x.sup.2] where [DELTA]t is the time step and [DELTA]x the spatial grid size, is critical for the numerical stability of the explicit finite difference schemes for this class of PDE, 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 PDE (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.

3 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 (DTM) 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.sup.k]}.sup.+[infinity].sub.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 DTM 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

[phi](t, k) = [d.sup.k]x (t)/d[t.sup.k], [for all]t [member of] D.

For t = [t.sub.i], then [phi](t, k) = [phi]([t.sub.i], k), where k belongs to a set of non-negative integers, denoted as the K domain. Therefore, (2) can be rewritten as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3)

where X(k) is called the spectrum of x(t) at t = [t.sub.i].

Definition 3.2. Suppose that x(t) is analytic in the time domain D, then it can be represented as

x(t) = [[infinity].summation over (k=0)] [(t - [t.sub.i]).sup.k]/k! X (k). (4)

Thus, the equation (4) represents the inverse transformation of X(k). Definition 3.3. If X(k) is defined as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

where k [member of] [Z.sup.+] [union] {0}, then the function x(t) can be described as

x(t) = 1/q(t) [[infinity].summation over (k=0)] [(t - [t.sub.i]).sup.k]/k! X(k)/M(k), (6)

where M(k) [not equal to] 0 and q(t) [not equal to] 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.sub.0], [t.sub.1], ..., [t.sub.n] = H} such that [t.sub.i] < [t.sub.i+1] and = [t.sub.i+1] - [t.sub.i] for i = 0, ..., n - 1.

Let M(k) = [H.sup.k.sub.i]/k!, q(t) = 1 and x(t) be a analytic function in [0, H], It then defines the differential transformation as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (7)

and its differential inverse transformation of X(k) is defined as follow

x(t) = [[infinity].summation over (k=0)] [(t/[H.sub.i]).sup.k] X(k), for t [member of] [[t.sub.i], [t.sub.i+1]]. (8)

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 DTM, 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.,

x(t) = 1/q(t) [n.summation over (k=0)] [(t - [t.sub.i]).sup.k]/k! X(k)/M(k) + [R.sub.n+1] = [n.summation over (k=0)] [(t/H).sup.k]X (k) + [R.sub.n+1], (9)

where

[R.sub.n+1] = [[infinity].summation over (k=n+1)] [(t/H).sup.k] X(k), and [R.sub.n+1] [right arrow] 0 as n [right arrow] [infinity].

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.sub.i], where H = [[summation].sup.n.sub.i=1] [H.sub.i]. It is important to remark that, [H.sub.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 [less than or equal to] t [less than or equal to] [t.sub.1], [t.sub.0] = 0), the one dimensional differential transformation is given by

x(t) = [n.summation over (k=0)] [(t/[H.sub.0]).sup.k] [X.sub.0] (k), where [X.sub.0](0) = [x.sub.0](0). (10)

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.sub.i] then one gets, step by step, an approximate solution in the whole domain [0, H], each initial value x([H.sub.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.sub.0] can be solved entirely in the first subdomain. The end point of function x(t) in the first subdomain is [x.sub.1], and the value of t is [H.sub.0]. Thus, [x.sub.1](t) is obtained by the differential transformation method as

[x.sub.1]([H.sub.0]) = x([H.sub.0]) = [n.summation over (k=0)] [X.sub.0](k). (11)

Since that [x.sub.1]([H.sub.0]) represents the initial condition in the second subdomain, then [X.sub.1](0) = [x.sub.1]([H.sub.0]). And so the function x(t) can be expressed in the second sub-domain as

[x.sub.2]([H.sub.1]) = x([H.sub.1]) = [n.summation over (k=0)] [X.sub.1](k). (12)

In general, in each i - 1 subdomain one gets,

[x.sub.i]([H.sub.i]) = [x.sub.i-1]([H.sub.i-1]) + [n.summation over (k=1)] [X.sub.i-1] (0) + [n.summation over (k=1)] [X.sub.i-1] (k), i = 1,2, ..., n. (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 DTM and finite differences for reaction-diffusion PDE

Our goal in this section is to construct a hybrid scheme for Eq. (1i) combining the DTM and finite differences as follows: we begin taking the differential transforms of both sides of the governing equations from the 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (14i)

U(0, k) = [G.sub.1](k), k [greater than or equal to] 0, (14ii)

U(1, k) = [G.sub.2](k), k [greater than or equal to] 0, (14iii)

U(x, 0) = g(x), 0 [less than or equal to] x [less than or equal to] 1, (14iv)

where U(x, k), [G.sub.1](k), [G.sub.2](k) are the differential transform of u(x, t), [g.sub.1](t), and [g.sub.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 [DELTA]x of the interval 0 [less than or equal to] x [greater than or equal to] 1. The x coordinates of the grid points are given by [x.sub.j] = j[DELTA]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

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (15i)

[U.sub.0](k) = [G.sub.1](k), k [greater than or equal to] 0, (15ii)

[U.sub.M](k) = [G.sub.2](k), k [greater than or equal to] 0, (15iii)

[U.sub.j](0) = [g.sub.j], For j = 0, ..., M, (15iv)

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.,

[u.sub.j](t) = [n.summation over (k=0)] [(t/[H.sub.i]).sup.k] [U.sub.j] (k), 0 [less than or equal to] t [less than or equal to] [H.sub.i], (16)

and then the solution is given by:

u(x, t) = [M-1.summation over (j=1)] [u.sub.j](t), for 0 < x < 1. (17)

5 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 PDE. 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 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 DTM is faster than the multistage Adomian method, but the Runke-Kutta methods require less computation time in comparison with DTM and multistage Adomian [12].

Example 5.1. We consider the following deterministic problem

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where u(x, t) is unknowns in (0,1) x (0,1}, and the terms f (x, t) = [e.sup.x+t] (1 - (1 + 0.1xt) - 0.1t) and p(x, t) = 1 + 0.1xt are given. The exact solution is given by u(x, t) = [e.sup.x+t], [30].

We also consider the random version of Example (5.1), by introducing the following perturbations:

u(x, 0, [gamma]) = [e.sup.x] + [gamma], (18i)

u(0, t, [gamma]) = [e.sup.t] + [gamma], u(1, t, [gamma]) = [e.sup.1+t] + [gamma], (18ii)

p(x, t, [gamma]) = 1 + 0.1xt + [gamma], and (18iii)

f (x, t, [gamma]) = - 0.1[e.sup.x](1 + x)t[e.sup.t] + [gamma], (18iv)

where the unknown u(x, t, [gamma]) as well as u(x, 0, [gamma]), u(0, t, [gamma]), u(1, t, [gamma]), p(x, t, [gamma]), f (x, t, [gamma]), are stochastic processes (s.p.) defined on a common probability space ([OMEGA], F, P) and [gamma] is the random variable (r.v.), which consists of three components:

1. The sample space [OMEGA] which is a non-empty set that collects all the elementary events or states [epsilon].

2. A collection of subsets of [OMEGA] called [sigma]-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 [gamma] is a random variable (r.v.) defined on the original sample space [OMEGA] with probability function.

A solution of Example (5.1) with randomness means that for each [gamma] [member of] [OMEGA], u(x, t, [gamma]) 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:

[U.sub.j] (k, [gamma]) =[e.sup.j[DELTA]x] + [delta](k)[gamma], (19i)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (19ii)

[P.sub.j[+ or -]1/2] (k, [gamma]) = [delta](k) + 0.1(j [+ or -] 1/2) [DELTA]x[delta](k - 1) + [delta](k)[gamma], and (19iii)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (19iv)

where i denotes the i-th split domain.

5.1 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 [DELTA]x = 0.05 is used. The accuracy of the DTM 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 [DELTA]x and a fixed time step size [DELTA]t = 0.0001. The accuracy of the DTM 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 PDE 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 PDE. 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.

[FIGURE 1 OMITTED]

5.2 Randomness on the initial condition

Graphical results of Monte Carlo simulations using 50 realizations for the random version of Example (5.1) presented in and assuming a uniform probabilistic distribution for the initial conditions are shown on the left hand side of Figure 2. It can be observed that the expected solution (i.e. average of the ensemble of the realizations) does not agree with the solution of the deterministic version Example (5.1). However, it can be seen on the right hand side that when the number of realizations is increased to 200 the expected solution agrees very well with the solution of the deterministic version of Example (5.1).

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.

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

5.3 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).

[FIGURE 4 OMITTED]

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 [gamma]. This is due to the fact that perturbing the boundary conditions leads to an exact solution u(x, t) = [e.sup.x+t] + [gamma].

5.4 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.

[FIGURE 5 OMITTED]

5.5 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.

[FIGURE 6 OMITTED]

6 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 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 DTM 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 PDE 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 PDE has been selected due to the fact that reaction-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 some probability distributions functions are assumed for the diffusion coefficient, source term, boundary and initial conditions.

The effect of introducing randomness in the diffusion PDE is justified by the fact that the diffusion coefficient, source term, boundary and initial conditions have some degree of uncertainty. Therefore, the random diffusion PDE 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 PDE. Future works can be developed for more complex cases such a nonlinear system with two interacting scalar fields.

Acknowledgements

The authors thank the anonymous reviewers for their helpful suggestions and remarks. The first author has been partially supported by CDCHTAULA project grant I-1289-11-05-A.

References

[1] F. Boano, R. Revelli, and L. Ridolfi, "Stochastic modelling of DO and BOD components in a stream with random inputs," Advances in Water Resources, vol. 29, no. 9, pp. 1341-1350, 2006. 51, 52

[2] B. Chen-Charpentier, B. Jensen, and P. Colberg, "Random Coefficient Differential Models of Growth of Anaerobic Photosynthetic Bacteria," ETNA, vol. 34, pp. 44-58, 2009. 51

[3] B. Oksendal, Stochastic Differential Equations. Springer, New York, 1995. 51

[4] T. Soong, Probabilistic Modeling and Analysis in Science and Engineering. Wiley, New York, 1992. 51, 62

[5] B. M. Chen-Charpentier, J.-C. Cortes, J.-V. Romero, and M.-D. Rosello, "Some recommendations for applying gPC (generalized polynomial chaos) to modeling: An analysis through the Airy random differential equation," Applied Mathematics and Computation, vol. 219, no. 9, pp. 4208-4218, 2013. 51

[6] G. Gonzalez-Parra, B. Chen-Charpentier, and A. J. Arenas, "Polynomial chaos for random fractional order differential equations," Applied Mathematics and Computation, vol. 226, pp. 123-130, 2014. 51

[7] E. Vanden Eijnden, "Studying random differential equations as a tool for turbulent diffusion," Phys. Rev. E, vol. 58, no. 5, pp. R5229-R5232, Nov 1998. [Online]. Available: http://link.aps.org/doi/10.1103/PhysRevE.58.R5229 51

[8] B. Kegan and R. West, "Modeling the simple epidemic with deterministic differential equations and random initial conditions," Math. Biosc, vol. 195, pp. 197-193, 2005. 51

[9] H. Kim, Y. Kim, and D. Yoon, "Dependence of polynomial chaos on random types of forces of KdV equations," Applied Mathematical Modelling, vol. 36, no. 7, pp. 3080-3093, 2012. 51

[10] J. Wu, Y. Zhang, L. Chen, and Z. Luo, "A Chebyshev interval method for nonlinear dynamic systems under uncertainty," Applied Mathematical Modelling, vol. 37, no. 6, pp. 4578-4591, 2013. 51

[11] S. Bhatnagar and Karmeshu, "Monte-Carlo estimation of time-dependent statistical characteristics of random dynamical systems," Applied Mathematical Modelling, vol. 35, no. 6, pp. 3063-3079, 2011. 51, 52

[12] G. Gonzalez-Parra, L. Acedo, and A. Arenas, "Accuracy of analytical-numerical solutions of the Michaelis-Menten equation," Computational & Applied Mathematics, vol. 30, no. 2, pp. 445-461, 2011. 51, 61

[13] L. Villafuerte and B. Chen-Charpentier, "A random differential transform method: Theory and applications," Applied Mathematics Letters, vol. 25, no. 10, pp. 1490-1494, 2012. 51

[14] F. Morrison, The Art of Modeling Dynamic Systems. John Wiley, 1991. 52

[15] Y. Chen, J. Liu, and G. Meng, "Incremental harmonic balance method for nonlinear flutter of an airfoil with uncertain-but-bounded parameters," Applied Mathematical Modelling, vol. 36, no. 2, pp. 657-667, 2012. 52

[16] V. Mallet and B. Sportisse, "Air quality modeling: From deterministic to stochastic approaches," Comput. Math. Appl, vol. 55, no. 10, pp. 2329-2337, 5 2008. 52

[17] S. D. Brown, R. Ratcliff, and P. L. Smith, "Evaluating methods for approximating stochastic differential equations," Journal of Mathematical Psychology, vol. 50, no. 4, pp. 402-410, 8 2006. 52

[18] S. Wu, "The Euler scheme for random impulsive differential equations," Applied Mathematics and Computation, vol. 191, no. 1, pp. 164-175, 2007. 52

[19] S. R. Hanna, J. C. Chang, and M. E. Fernau, "Monte Carlo estimates of uncertainties in predictions by a photochemical grid model (uam-iv) due to uncertainties in input variables," Atmospheric Environment, vol. 32, no. 21, pp. 3619-3628, 1998. 52, 62

[20] K. M. Hanson, "A framework for assessing uncertainties in simulation predictions," Physica D: Nonlinear Phenomena, vol. 133, no. 1-4, pp. 179-188, 1999. 52

[21] G. Pukhov, Differential Transformations of Functions and Equations. Naukova Dumka (in Russian), 1980. 54, 57

[22] J. Zhou, Differential Transformation and its Applications for Electrical Circuits. Huazhong University Press, Wuhan (in Chinese), 1986. 55, 57

[23] J. Biazar and M. Eslami, "Differential transform method for quadratic Riccati differential equation," International Journal of Nonlinear Science, vol. 9, no. 4, pp. 444-447, 2010. 55

[24] A. J. Arenas, G. Gonzalez-Parra, and B. M. Chen-Charpentier, "Dynamical analysis of the transmission of seasonal diseases using the differential transformation method," Mathematical and Computer Modelling, vol. 50, no. 5-6, pp. 765-776, 2009. [Online]. Available: http://dx.doi.org/10.10167j.mcm.2009.05.005 55, 57

[25] I. A.-H. Hassan, "Application to differential transformation method for solving systems of differential equations," Applied Mathematical Modelling, vol. 32, pp. 2552-2559, 2008. 55

[26] M. Jang and C. Chen, "Analysis of the response of a strongly nonlinear damped system using a differential transformation technique," Applied Mathematics and Computation, vol. 88, pp. 137-151, 1997. 55

[27] V. S. Erturk, G. Zaman, and S. Momani, "A numeric-analytic method for approximating a giving up smoking model containing fractional derivatives," Computers & Mathematics with Applications, vol. 64, no. 10, pp. 3065-3074, 2012. 55

[28] C. Bervillier, "Status of the differential transformation method," Applied Mathematics and Computation, vol. 218, no. 20, pp. 10 158-10 170, 2012. 55, 57

[29] I. Hwang, J. Li, and D. Du, "A numerical algorithm for optimal control of a class of hybrid systems: differential transformation based approach," International Journal of Control, vol. 81, no. 2, pp. 277-293, 2008. [Online]. Available: http://dx.doi.org/10.1080/00207170701556880 55

[30] C. E. Mejia and D. A. Murio, "Numerical identification of diffusivity coefficient and initial condition by discrete mollification," Comput. Math. Appl, vol. 12, pp. 35-50, 1995. 61

[31] C.-K. Chen and S.-P. Ju, "Application of differential transformation to transient advective-dispersive transport equation," Applied Mathematics and Computation, vol. 155, no. 1, pp. 25-38, 2004. 63

Gilberto Gonzalez-Parra (1), Abraham J. Arenas (2) and Miladys Cogollo (3)

(1) Universidad de Los Andes, Merida, Venezuela, gcarlos@ula.ve.

(2) Universidad de Cordoba Monteria, Colombia, aarenas@correo.unicordoba.edu.co.

(3) Universidad EAFIT, Medellin, Colombia, mcogollo@eafit.edu.co.

Received: 09-02-2015 | Accepted: 18-06-2015 | Online: 31-07-2015

MSC: 35K10, 65C05, 65M99, 35R60

doi: 10.17230/ingciencia.11.22.3
```Table 1: Numerical absolute errors using the hybrid method with 3-
terms, diffe-rent time step sizes and a space step size [DELTA]x =
0.05 for the deterministic version of Example (5.1) at x = 0.5

Time   h = 0.001   h = 0.0005

0.1    0.39E--4    0.52E--5
0.2    0.58E--4    0.70E--5
0.5    0.82E--4    0.61E--5
1.0    0.12E--3    0.31E--5
2.0    0.24E--3    0.76E--4

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

Time   [DELTA]x = 0.05   [DELTA]x = 0.04   [DELTA]x = 0.033

0.1       0.39E--4          0.11E--4           0.59E--5
0.2       0.58E--4          0.17E--4           0.93E--5
0.5       0.82E--4          0.29E--4           0.16E--4
1.0       0.12E--3          0.56E--4           0.31E--4
2.0       0.24E--3          1.80E-04           0.11E--3
```