ARTÍCULO ORIGINAL
DOI: 10.17230/ingciencia.11.22.3
Analytical-Numerical Solution of a Parabolic Diffusion Equation Under Uncertainty Conditions Using DTM with Monte Carlo Simulations
Solucion numérico-analítica de una ecuación de difusión bajo condiciones de incertidumbre utilizando DTM y Monte Carlo
Gilberto González-Parra1, Abraham J. Arenas2 and Miladys Cogollo 3
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 .
Received: 09-02-2015 | Accepted: 18-06-2015 | Online: 31-07-2015
MSC: 35K10, 65C05, 65M99, 35R60
Abstract
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 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.
Key words: random linear diffusion models; uncertainty conditions; finite difference schemes; differential transformation method; analytical-numerical solution
Resumen
Un método numérico para resolver una ecuación parabólica general aleatoria lineal donde el coeficiente de difusión, el término fuente, las condiciones de contorno e iniciales incluyen la incertidumbre, se ha desarrollado. Ecuaciones de difusión surgen en muchos campos de la ciencia y la ingeniería, y en muchos casos, existen la incertidumbres debido a los datos que no se pueden saber, o debido a errores en las mediciones y la variabilidad intrínseca. Para modelar estas incertidumbres los parámetros correspondientes, coeficiente de difusión, término fuente, condiciones de contorno e iniciales, se suponen que son variables aleatorias con determinadas distribuciones de probabilidad. Basándose en los resultados numéricos, se obtienen los intervalos de confianza y valores medios esperados para la solución. Además, se obtienen con las soluciones numéricas-analíticas del método híbrido propuesto.
Palabras clave: modelos lineales de difusión aleatorios; condiciones de incertidumbre; Esquemas de diferencias finitas; método de transformación diferencial; solución analítica-numérica
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:
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:
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 , 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 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 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
For t=ti, then φ(t,k)=φ(ti,k), where k belongs to a set of non-negative integers, denoted as the K domain. Therefore, (2) can be rewritten as
where X(k) is called the spectrum of x(t) at t=ti.
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).
Definition 3.3 If X(k) is defined as
where k
+ ∪{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=t0,t1,... ,tn=H} such that
ti < ti+1 and Hi=ti+1−ti for i=0,... ,n−1.
Let M(k)= , q(t)=1 and x(t) be a analytic function
in [0,H]. It then defines the differential transformation as
and its differential inverse transformation of X(k) is defined as follow
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.,
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 Hi,
where H=. It is important to remark that,
Hi 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 ≤ t1, t0=0), the one dimensional
differential transformation is given by
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 ti then one gets, step by step, an approximate solution in the whole domain [0,H], each initial value x(Hi 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 X0 can be solved entirely in the first subdomain. The end point of function x(t) in the first subdomain is x1, and the value of t is H0. Thus, x1(t) is obtained by the differential transformation method as
Since that x1(H0) represents the initial condition in the second subdomain, then X1(0)=x1(H0). And so the function x(t) can be expressed in the second sub-domain as
In general, in each i−1 subdomain one gets,
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:
where U(x,k), G1(k), G2(k) are the differential transform of u(x,t), g1(t), and g2(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 xj=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:
where bold uppercase letters are for the transformed functions. Thus, from a process of inverse differential transformation, the solutions on each sub-domain can be obtained using n+1 terms for the power series, i.e.,
and then the solution is given by:
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
where u(x,t) is unknowns in (0,1)×(0,1], and the terms f(x,t)=ex+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)=ex+t, [30].
We also consider the random version of Example (5.1), by introducing the following perturbations:
where the unknown u(x,t,γ) as well as u(x,0,γ), u(0,t,γ), u(1,t,γ), p(x,t,γ),
f(x,t,γ), are stochastic processes (s.p.) defined on a common probability space (Ω,,P) and γ is the random variable (r.v.),
which consists of three components:
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.
5.1 Numerical comparisons in accurac
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 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 ∆x and a fixed time step size ∆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.
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.
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).
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)=ex+t+γ.
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.
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.
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 CDCHTA-ULA 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.
[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.
[3] B. Oksendal, Stochastic Differential Equations. Springer, New York, 1995.
[4] T. Soong, Probabilistic Modeling and Analysis in Science and Engineering. Wiley, New York, 1992.
[5] B. M. Chen-Charpentier, J.-C. Cortés, J.-V. Romero, and M.-D. Roselló, ''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.
[6] G. González-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.
[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
[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.
[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.
[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.
[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.
[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.
[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.
[14] F. Morrison, The Art of Modeling Dynamic Systems. John Wiley, 1991.
[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.
[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.
[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.
[18] S. Wu, ''The Euler scheme for random impulsive differential equations,'' Applied Mathematics and Computation, vol. 191, no. 1, pp. 164-175, 2007.
[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.
[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.
[21] G. Pukhov, Differential Transformations of Functions and Equations. Naukova Dumka (in Russian), 1980.
[22] J. Zhou, Differential Transformation and its Applications for Electrical Circuits. Huazhong University Press, Wuhan (in Chinese), 1986.
[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.
[24] A. J. Arenas, G. González-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.1016/j.mcm.2009.05.005
[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.
[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.
[27] V. S. Ertürk, 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.
[28] C. Bervillier, ''Status of the differential transformation method,'' Applied Mathematics and Computation, vol. 218, no. 20, pp. 10 158 - 10 170, 2012.
[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
[30] C. E. Mejía and D. A. Murio, ''Numerical identification of diffusivity coefficient and initial condition by discrete mollification,'' Comput. Math. Appl, vol. 12, pp. 35-50, 1995.
[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.