Abstract
The consideration of uncertainty in differential equations leads to the emergent area of random differential equations. Under this approach, inputs become random variables and/or stochastic processes. Often one assumes that inputs are independent, a hypothesis that simplifies the mathematical treatment although it could not be met in applications. In this paper, we analyse, through the Airy equation, the influence of statistical dependence of inputs on the output, computing its expectation and standard deviation by Fröbenius and Polynomial Chaos methods. The results are compared with Monte Carlo sampling. The analysis is conducted by the Airy equation since, as in the deterministic scenario its solutions are highly oscillatory, it is expected that differences will be better highlighted. To illustrate our study, and motivated by the ubiquity of Gaussian random variables in numerous practical problems, we assume that inputs follow a multivariate Gaussian distribution throughout the paper. The application of Fröbenius method to solve Airy equation is based on an extension of the method to the case where inputs are dependent. The numerical results show that the existence of statistical dependence among the inputs and its magnitude entails changes on the variability of the output.
1. Introduction and Motivation
Deterministic differential equations (ddes) have demonstrated to be powerful tools for modelling numerous problems appearing in different areas including physics, chemistry, economy, engineering, and epidemiology. Their practical application requires knowing their inputs (coefficients, forcing terms, initial/boundary conditions, etc.). This task can only be done after accurate measurements that usually contain uncertainty due to measuring errors or the inherent complexity or ignorance of the phenomena under study. This approach leads us to consider the inputs of such models as random variables (rvs) or stochastic processes (sps) rather than deterministic constants or functions, respectively. Differential equations containing in their formulation randomness are usually referred to as random differential equations (rdes) or stochastic differential equations (sdes) depending on the kind of uncertainty therein. When randomness is just considered through the white noise process (i.e., the generalized derivative of the Wiener process), they are usually called sdes Then, Itô calculus is required in order to conduct the study. Otherwise, the term rde is used (see, [1, page 66], [2]).
Rdes constitute a natural extension of ddes. Generalized Polynomial Chaos (usually denoted by gPC) and Monte Carlo sampling (MCs) are probably the most popular techniques to deal with rdes (see for instance [3, 4], resp.). In addition to these approaches, the extension of some deterministic techniques to the random scenario, based on the so-called -stochastic calculus, also constitutes useful tools to solve rdes (see [5, 6] and the references therein). In particular, a random power series Fröbenius method has been recently proposed by some of the authors to study some significant rdes by assuming that random inputs are independent [7, 8]. Although independence assumption simplifies the mathematical treatment of the models, it should not be assumed in many practical situations. Apart from few contributions such as [9, 10] where the authors study the dependent scenario by taking advantage of gPC approach, most methods developed to study rdes rely on independence of random inputs. In particular, to the best of our knowledge, applications of the random Fröbenius method considering dependent rvs have not been studied yet. As a consequence, the study of rdes with dependent inputs is currently an active research area, mainly stimulated by the necessity of providing more realistic approaches and more accurate answers in mathematical modelling.
In this paper, we present a comparative study about the capability of previous approaches to deal with rdes containing dependent random inputs. To conduct the study, we will consider the Airy random differential equation: where are assumed to be Gaussian dependent rvs on a probability space . We point out that the Airy equation has been selected since, as it is well-known, in the deterministic scenario its solutions are highly oscillatory [11]; therefore, it is expected that differences among gPC, Fröbenius method, and MCs will be better highlighted.
Specifically, we will compare, by means of several illustrative examples, the quality of the numerical approximations provided by the three approaches to compute the average and standard deviation of the solution sp to the initial value problem (ivp) (1). These examples will allow us to elucidate, through the random Airy differential equation, whether the statistical independence between the random inputs (initial conditions and coefficients), usually assumed in many applications, has a significant influence on the output.
The paper is organized as follows. Section 2 is divided in two parts. The first one is devoted to construct a mean square convergent random power series solution to the ivp (1) using the random Fröbenius method including approximations of the average and the variance of the solution sp In the second part, we summarize the main features of the gPC method to study rdes, and we apply gPC to the particular case where inputs are Gaussian dependent rvs. In Section 3 we will present several illustrative examples. The aim of these examples is twofold. First, to highlight the similarities and differences between the three approaches in dealing with both dependent/independent (Gaussian) rvs; second, to reveal, through the Airy rde, the importance of setting appropriately the statistical dependence of random inputs in dealing with mathematical models. Conclusions are drawn in the last section.
2. Development
In this section, we present the main results required to construct using Fröbenius and gPC methods the approximations of the mean and standard deviation of the solution sp to the ivp (1) when inputs ,, are assumed to be dependent rvs. We point out that the description of Fröbenius method in this scenario is presented briefly deliberately since it follows in broad outline that of independent case, which has been developed by some of the authors previously [7]. Foundations and further details about gPC method can be found in [3], for instance.
2.1. Tackling Random Dependence Using Fröbenius Method
Random Fröbenius method consists of constructing a power series solution to the ivp (1), say, which is mean square convergent on a certain -domain. The rigorous construction of such a mean square convergent random infinite series requires -calculus with , where convergence is defined in the -norm: [5, 7]. Convergence with /-norm is usually referred to as mean square (ms)/fourth (mf) convergence. By the Schwarz inequality, it is straightforward to proof that mean fourth convergence implies mean square convergence. Based on the ideas developed in [7], now we assume that the following joint absolute moments with respect to the origin increase at most exponentially; that is, there is a nonnegative integer and positive constants and , such that In the case of ivp (1), the random power series solution is given by [7] Under hypothesis (2) mf convergence (and hence ms convergence) of the first series appearing in (3) follows straightforwardly: where . As a consequence, for we have obtained as a majorant series that is convergent for all as it can be directly checked by D’Alembert test: The mf convergence for the second series in (3) follows analogously.
Remark 1. By assuming that there are positive constants and such that for every and for , m.s. convergence of series (3) can also be established analogously as we have shown previously. In fact, this follows immediately from the Schwarz inequality where , being < , . However, notice that this condition is stronger than (2).
Taking advantage of m.s. convergence of series appearing in right-hand side of (3) together with the following property ([2, page 88]): we can obtain approximations for the average, , and variance, , (or equivalently standard deviation) by truncating the random power series of given by (3). For the approximation of the average one gets In order to obtain approximations of the standard deviation of , we will take into account the well-known representation of the variance in terms of the two first moments, . Therefore, it is enough to approximate the second moment:
Remark 2. In order to legitimate the use of the previous approximations to the average and the standard deviation, condition (2) must be checked in practice. However, there is a lack of explicit formulae for the absolute moments with respect to the origin of some rvs. This aims us to look for a general approach to deal with a wide range of random inputs taking advantage of the so-called censuring method (see [12, chapter V]). Let us assume that rvs , , satisfy Then where denotes the joint probability density function (p.d.f.) of rvs , and , , being , . Indeed, in the case that , one gets Notice that in the last step the double integral is just 1, since is a pdf The other cases can be analyzed analogously. Substituting the integral by a sum in (12), previous reasoning remains true when and/or , are discrete rvs. As a consequence, important rvs such as binomial, hypergeometric, uniform or beta satisfy condition (2), which is related to joint absolute moments of , , . It is worthwhile to point out that there are significant rvs that do not satisfy condition (2) such as the exponential rv, for instance. In fact, taking , in (2), if , , then . As a consequence, in this case condition (2) does not fulfill. Although other unbounded rvs can also verify condition (2), we do not need to check it each time, since if we censure its codomain suitably, we are legitimated to compute approximations to the mean and standard deviation according to formulae (9)-(10), respectively. The larger the censured interval, the better the approximations. However, in practice, intervals relatively short provide very good approximations. For instance, as an illustrative example notice that the truncated interval contains the of the probability mass of a Gaussian rv with mean and standard deviation .
2.2. Tackling Random Dependence by Generalized Polynomial Chaos Method
As it has been underlined in Section 1, gPC constitutes a powerful method to deal with randomness in differential equations, say where denotes a differential operator; is the solution sp to be determined and is a forcing term. Notice that in the rde (14) uncertainty is represented by and it just enters through its coefficients and forcing term, although in practice it could also be considered via initial and/or boundary conditions. For the sake of clarity, in the following each scalar random input will be denoted by .
gPC permits to represent spectrally each in the random dimension, and the solution sp, , in . These representations are given by infinite random series defined in terms of certain orthogonal polynomial expansions which depend on a number of rvs , , This set constitutes a complete orthogonal basis in with the inner product where is the Kronecker delta and denotes the ensemble average defined as follows: being the joint pdf of and its support.
The choice of the trial basis is crucial in dealing with rdes. In [3], authors provide a comprehensive way to choose the trial basis according to the statistical distribution of the random input in order to achieve optimal convergence in (15). For instance, if rv follows a binomial, negative binomial, hypergeometric, Poisson, beta, or gamma distribution, then should be taken as Krawtchouk, Meixner, Hahn, Charlier, Jacobi, Laguerre orthogonal polynomials belonging to the Wiener-Askey scheme, respectively. In the significant case that is a Gaussian rv, Hermite polynomials are required. This particular case is referred to as Polynomial Chaos rather than gPC. Throughout this paper only PC will be used. The key connection to do an adequate selection of the trial basis lies in the close relationship between the pdf of some standard rvs and the weight function that defines the inner product (17) with respect to which some classical polynomials are orthogonal.
In order to keep the computations affordable in dealing with rdes, each random model parameter as well as the solution sp is represented by truncated series of the form (15), where the number of components of random vector also needs to be truncated at a number called the order of chaos, . The truncation order is made so that all expansion polynomials up to a certain maximum degree, , are included. This entails the following relationship between the number of terms in the series expansions, the maximum degree , and the order of chaos : In this context, solving the rde (14) consists of computing coefficients appearing in (18) which also allows to compute approximations of the expectation and the standard deviation to the solution sp as follows: To achieve this goal, first the expansion of given by (18) is substituted into the rde (14). Second, a Galerkin projection is done by multiplying the rde by every polynomial of the expansion basis , and then, the ensemble average is taken. This leads to that corresponds to a deterministic system of coupled differential equations whose unknowns are the node functions . These unknowns can be computed by standard numerical techniques such as Runge-Kutta scheme.
Most of the contributions based on gPC assume that rvs , are independent which facilitates the study. The case in which random parameters are assumed to be dependent is currently a topic under study. In [9, 10], authors present methods based on gPC to tackle dependence in differential equations. Both contributions provide general techniques that can be applied whenever the joint pdf of the random inputs is known. However, in practice the availability of this joint pdf can be very difficult even impossible. In the particular case where the inputs are dependent Gaussian rvs, an alternative method can be applied to conduct the corresponding study taking advantage that uncorrelation and independence are equivalent notions for Gaussian rvs together with Cholesky matrix decomposition. To exhibit how the method is going to be applied in our case, let us remember the following basic result.
Proposition 3. Let be a Gaussian vector with mean and variance-covariance matrix : (the symbol ⊤ denotes the usual matrix transpose operator). For each deterministic vector and deterministic matrix , the random vector defined by the linear transformation: follows the Gaussian distribution: .
On the other hand, let be a Gaussian random vector. As is a variance-covariance matrix, it is Hermitian and positive definite. Hence, there is a matrix, say , such that . For instance, Cholesky decomposition provides a well-known procedure to compute matrix [13].
Keeping the notation of the previous context, we apply Proposition 3 to the particular case where (so, and , being the identity matrix of size ); and . Then As , then , are Gaussian and uncorrelated r.v.’s and, therefore, independent. As a consequence, expression (22) provides a direct way to represent a Gaussian random vector with components statistically dependent by means of a linear transformation of a Gaussian vector whose components are independent: .
Now we detail how the previous development can be applied to transform the ivp (1), where random inputs ,, are assumed to be Gaussian dependent into another one with Gaussian independent random inputs. Let be the multivariate Gaussian distribution of the random data, and let us denote by the Cholesky decomposition of variance-covariance matrix of . According to (22), we define the linear transformation: where denotes the mean of vector and , being independent and identically distributed standard Gaussian rvs, that is, , . By (23), ivp (1) can be recast as follows: where random inputs , , are independent standard Gaussian rvs. This allows us to compute approximations of the expectation and standard deviation functions by PC according to (20).
3. Examples
In this section we will present several illustrative examples based on ivp (1) in order to compare the approximations provided by Polynomial Chaos, Fröbenius, and Monte Carlo simulation. The comparison is performed by computing the average and standard deviation functions of the solution of ivp (1). As we pointed out in the Introduction section, the study is conducted through Airy differential equation since, as in the deterministic case its solutions are highly oscillatory, it is expected that differences among the three previous approaches will be better highlighted in the random framework. Computations have been carried out with Mathematica package [14]. In particular, the coupled systems of differential equations obtained after applying gPC in each example are numerically solved with this software.
The examples have been designed to explore both the marginal influence of randomness on the output when it is assumed that only some inputs of ivp (1) are rvs (see Examples 1 and 2, where and are assumed to follow a bivariate Gaussian rvs, resp.), and the cases presented in Examples 3 and 4 where all the random inputs are assumed to follow a multivariate Gaussian distribution. In the two first examples, we also investigate the influence of the numerical value of the correlation coefficient of the two-dimensional random input changes. Examples 3 and 4 seek to illustrate the different qualitative behaviour of the solution sp of ivp (1) depending on rv .
Example 1. This first example has been devised to investigate whether statistical dependence between the initial conditions entails a substantial change in the output with respect to independence assumption. In addition, it permits to highlight some significant advantages of Fröbenius and PC methods in comparison with MCs. Let us consider the ivp (1) where and the initial conditions are assumed to be dependent on a Gaussian r.v.’s: , where
In this case, we can directly monitor the influence of the statistical dependence between and , measured through its correlation coefficient , in the computations of the average and the standard deviation when these moments are calculated using anyone of the three methods. For the Fröbenius method, taking into account that and expression (10), we observe that dependence only contributes through the term: . In case of PC method, monitoring dependence is made directly over the resulting ivp:
since some coefficients , , depend on . Specifically, are the entries of the Cholesky decomposition of variance-covariance matrix . In this case one gets
Notice that in (26), and are independent standard Gaussian rvs: . With respect to MC simulation, it is obvious that dependence is monitored directly through the bivariate pdf which depends on :
The random samples needed to apply MC method have been sampled with the Mathematica instruction: “RandomVariate[MultinormalDistribution].”
Notice that Fröbenius and PC approaches turn out much more fruitful than MCs in this particular example. In the case of Fröbenius method, the series representation of the solution sp given by (3) permits not only to compute reliable approximations of the average and standard deviation but also to determine its full statistical distribution, which in general is a major challenge. Firstly, notice that under condition (2), we have proven that the infinite series (3) converges in the mean square sense for each ; therefore, it also converges in distribution for every . In accordance with (3), the truncated approximation of is given by
where we adopt the convention for . Note that
where denotes the confluent hypergeometric function of order and evaluated at the point .
Since the vector has a bivariate Gaussian distribution, follows a univariate Gaussian distribution whose average and variance are given by
respectively. By property (8), these expressions converge as to the exact average, , and variance , respectively. This determines completely the statistical distribution of for every .
On the other hand, PC method also provides, in this example, a useful series representation of the solution sp that permits to obtain its statistical distribution. It is straightforward to check that the PC series representation obtained from ivp (26) has the following finite linear form:
where , independent, and , and . Hence, following an analogous reasoning as in the previous case, we deduce the full distribution of . With respect to MC approach, it only provides a set of numerical values to the solution at some selected time instants from which only rough approximations of the statistical distribution of the solution can be given.
In Tables 1 and 2, we compare the approximations of the expectation and standard deviation at different time instants and correlation values: , respectively. Notice that these values correspond to average negative dependence; independence; average positive dependence; and strong positive dependence, respectively. We assume , , , . (), () and () denote the approximate expectation (and standard deviation) at time obtained by Fröbenius (see expressions (9) and (10), or equivalently, (31) and (32)) with truncation orden , PC (see expression (20)) of order and MCs using simulations, respectively. The values of and show in Tables 1 and 2 are those obtained when the numerical stabilization at six significant digits is achieved.
We observe the numerical results provided by Fröbenius and PC approaches match while MCs captures about three significant digits. From Tables 1 and 2, we realize that the correlation value between and does not influence the average, but it does decisively influuence the standard deviation of the solution. This indicates to us that it is crucial to know not only the existence of statistical dependence between initial conditions but also quantifying, as accurate as possible, its value by the correlation coefficient . Notice that these numerical conclusions agree with formulae (31)-(32).
Example 2. In the previous example, uncertainty entered in the equation just through initial conditions and . Although both rvs are dependent, we have showed that it does not influence the expectation of the solution but its standard deviation. Does this answer change in case that randomness is considered through coefficient and an initial condition? In order to answer this question let us consider the ivp (1) where are assumed to be dependent Gaussian rvs: , where
To perform the computations of the average () and standard deviation () shown in Figures 1 and 2, respectively, we have taken , , and . The initial condition is assumed to be deterministic: . Computations for Fröbenius and PC have been carried out until the numerical stabilization at six significant digits is achieved. This stabilization is got for Fröbenius method with for the average and for the standard deviation. The corresponding values for PC are and . These results for and do not depend on the correlation coefficient . Results obtained by MCs have been carried out with simulations. As a reference of the computational burden required by each one of the three methods, we indicate the CPU seconds to carry out the standard deviation plotted in Figure 2 in an Intel Core 7 with GHz: Fröbenius ( s), PC ( s), and MC ( s).
Although similar numerical differences between MCs and the other techniques (Fröbenius and PC) could be reported in a table as we did in the foregoing example, for both the average and the standard deviation, now we present the numerical approximations in Figures 1 and 2 without labelling each method since they are indistinguishable graphically. We underline that numerical results for Fröbenius and PC methods match with the six significant digits. While, MCs captures three of these digits.
In contrast to what happened in Example 1, we now observe that average changes, but slightly, when does (see Figure 1). These changes are greater when computing standard deviation (see Figure 2). Both conclusions agree with formulae (9)-(10). Again, these results indicate that independence between random parameters (in this case, between coefficient and the initial condition ) must be checked previously since the existence of statistical dependence influences significantly on the output.
Example 3. Now, we consider the case where all the data involved in ivp (1) are rvs. We will assume that the random vector follows a multivariate Gaussian distribution. The aim of this example is twofold: first, to confirm the conclusions drawn in the two previous examples and, second, to show and compare the capabilities of the three methods under analysis to tackling satisfactorily full randomness in the ivp (1). Computations have been performed taking as average: and, and as the variance-covariance matrices given by (35) that correspond to statistical independence and dependence, respectively,
Computations for Fröbenius and PC have been carried out until the numerical stabilization at six significant digits is achieved. This stabilization is got for Fröbenius method with for the average and for the standard deviation. The corresponding values for PC are and . As it happened in the foregoing example, the obtained results for and do not depend on the variance-covariance matrix. Numerical values obtained by MCs have been carried out with simulations. Again, as in the previous example the numerical results provided by Fröbenius and PC methods coincide with the six significant digits, while MCs only captures three of these digits.
Figures 3 and 4 show the results for the average and the standard deviation on the time interval , respectively. Analogous comments as we did in Example 2 can be done: one presents changes for both the expectation and standard deviation of the solution depending on the statistical dependence of the input, being these changes greater on the standard deviation.
Example 4. The aim of this example is to show that the qualitative behaviour of the solution sp of the random Airy differential (1) is different depending on the random input . To illustrate this fact, we will assume that the random vector follows a multivariate Gaussian distribution with average: and the same variance-covariance matrices as the ones considered in Example 3 (see expression (35)). Note that we are assuming that most part of the probability mass of rv is on the negative real line which implies a different behaviour of the solution with respect to previous examples. Figure 5 shows the mean of the solution in both cases, that is, when the random inputs are independent and dependent. Although both plots are quite similar, in Figure 6 we see in contrast to what happens in previous examples that now standard deviations as well as their difference increase over time.
(a)
(b)
4. Conclusions
The consideration of uncertainty in models based on differential equations leads to random differential equations. Over the last few decades, these types of random continuous models have demonstrated to be powerful tools in dealing with mathematical modelling. However, for simplicity, most of these contributions rely on the assumption that random inputs are statistically independent, a hypothesis that could not be met in many applications. The study of differential equations whose inputs are statistically dependent constitutes currently a topic under development.
In this paper we have studied two methods for solving differential equations whose parameters are assumed to be Gaussian dependent, namely, Fröbenius and Polynomial Chaos. The numerical results, for the average and standard deviation of the solution, provided by both methods have been compared with those computed by Monte Carlo simulations, which can be considered the most common approach to deal with random differential equations. The study has been performed through the Airy equation, which is expected to be an excellent test model to highlight differences among previous approaches, due to the highly oscillatory behaviour of its solutions in the deterministic case. The examples reveal that Fröbenius and Polynomial Chaos perform better than Monte Carlo simulations since both are more accurate. A major conclusion drawn from the study case performed through the examples is the significant influence of statistical dependence among the random inputs on the variability of the output. As a consequence, the usual hypothesis of statistical independence for the random parameters should be checked carefully in modelling. Furthermore, when dependence is assumed, its numerical value (measured, for example, by the correlation coefficient) must be determined as accurately as possible, since it has been shown that it also influences both the average and the variability of the output.
Notice that in order to conduct the study, we have had to extend the Fröbenius method presented by some of the authors in the previous contribution [7] to the case where random inputs are dependent. Although we have focused the study on the case in which random variables are just Gaussian, notice that Fröbenius approach does not depend on the statistical type of the involved random variables, while tackling statistical dependence with Polynomial Chaos has been carried out through a direct approach based on nice properties of Gaussian random variables. This approach has allowed us to transform the random initial value problem (1) into another having independent Gaussian random variables, which has facilitated the study.
Solving random differential equations mainly consists of computing the average and standard deviation of the solution stochastic process. A major challenge is to determine its statistical distribution. Example 1 shows, by means of a simple but still illustrative scenario, the potentiality of both, Fröbenius and gPC methods, to deal with this issue when other random differential equations appear in modelling. We think that the combined application of the novel theory of copulas [15] and previous methods constitutes a promising approach that will be considered in the forthcoming works.
Acknowledgments
This work has been partially supported by the Ministerio de Economía y Competitividad Grants MTM2009-08587 and DPI2010-20891-C02-01 and Universitat Politècnica de València Grant PAID06-11-2070.