- Email: [email protected]
A flexible numerical approach for quantification of epistemic uncertainty
Journal of Computational Physics 240 (2013) 211–224
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
A flexible numerical approach for quantification of epistemic uncertainty Xiaoxiao Chen a, Eun-Jae Park b,1, Dongbin Xiu a,c,⇑,2 a
Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA Department of Computational Science & Engineering, Yonsei University, Seoul 120-749, Republic of Korea c Department of Mathematics, and Scientific Computing and Imagining Institute University of Utah, Salt Lake City, UT 84112, USA b
a r t i c l e
i n f o
Article history: Received 2 December 2011 Received in revised form 15 January 2013 Accepted 19 January 2013 Available online 4 February 2013 Keywords: Uncertainty quantification Epistemic uncertainty Generalized polynomial chaos Encapsulation problem
a b s t r a c t In the field of uncertainty quantification (UQ), epistemic uncertainty often refers to the kind of uncertainty whose complete probabilistic description is not available, largely due to our lack of knowledge about the uncertainty. Quantification of the impacts of epistemic uncertainty is naturally difficult, because most of the existing stochastic tools rely on the specification of the probability distributions and thus do not readily apply to epistemic uncertainty. And there have been few studies and methods to deal with epistemic uncertainty. A recent work can be found in [J. Jakeman, M. Eldred, D. Xiu, Numerical approach for quantification of epistemic uncertainty, J. Comput. Phys. 229 (2010) 4648–4663], where a framework for numerical treatment of epistemic uncertainty was proposed. The method is based on solving an encapsulation problem, without using any probability information, in a hypercube that encapsulates the unknown epistemic probability space. If more probabilistic information about the epistemic variables is known a posteriori, the solution statistics can then be evaluated at post-process steps. In this paper, we present a new method, similar to that of Jakeman et al. but significantly extending its capabilities. Most notably, the new method (1) does not require the encapsulation problem to be in a bounded domain such as a hypercube; (2) does not require the solution of the encapsulation problem to converge point-wise. In the current formulation, the encapsulation problem could reside in an unbounded domain, and more importantly, its numerical approximation could be sought in Lp norm. These features thus make the new approach more flexible and amicable to practical implementation. Both the mathematical framework and numerical analysis are presented to demonstrate the effectiveness of the new approach. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Extensive research efforts have been devoted to uncertainty quantification (UQ), especially for large and complex systems. One of the key issues in UQ is to design efficient numerical algorithms to propagate uncertainty from system inputs to outputs. To this end, the most widely used approach adopts probabilistic framework, where the uncertainty inputs are modeled as random variables/processes. One condition naturally required by most probabilistic methods is the knowledge ⇑ Corresponding author. Tel.: +1 765 496 2846. E-mail addresses: [email protected] (E.-J. Park), [email protected]; [email protected] (D. Xiu). The work of this author was supported in part by the WCU program through the Korea Science and Engineering Foundation funded by the Ministry of Education, Science, and Technology R31-2008-000-10049-0. 2 The work of D. Xiu was supported by AFOSR, DOE, NNSA, and NSF; also by NRF R31-2008-000-10049-0 through Department of Computational Science & Engineering, Yonsei University, Seoul 120-749, Korea. Tel.: +1 202 642 4808. 1
0021-9991/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.01.018
212
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
of the probability distribution of the uncertain inputs. In the field of UQ, this kind of uncertainty is sometimes referred to as aleatory uncertainty — the kind of uncertainty whose complete probabilistic specification is available. In practice, however, it is often very difficult, if not impossible, to have sufficient information to fully prescribe the probability distribution of all the random inputs. And in such cases, one is facing epistemic uncertainty, the kind of uncertainty whose probability distribution is not known completely. Until recently, most uncertainty analysis has focused on aleatory uncertainty. Numerous stochastic methods have been developed to provide accurate and efficient simulations for this form of uncertainty. For example, the stochastic methods based on generalized polynomial chaos (gPC) [16], an extension of the traditional polynomial chaos [6], have demonstrated their efficiency in practice and are widely used for many problems. Their implementations typically follow either stochastic Galerkin (SG) [2,6,16] or stochastic collocation (SC) [1,4,12,15], often termed as intrusive methods or non-intrusive methods, respectively. For a detailed review on the methods, see [13,14]. Numerical study of epistemic uncertainty is more difficult because of the lack of sufficient probabilistic information. And one naturally can not readily adopt probabilistic approaches. Some of the existing approaches include evidence theory [8], possibility theory [3] and interval analysis [7,10]. These method have their own advantages, though most do not address efficient numerical implementations. More recent studies can be found in [5,11]. This paper is based on a more recent work of [9], where a general numerical approach was proposed for epistemic uncertainty analysis. The method consists of three components: (1) encapsulation of the epistemic uncertain inputs by a bounded domain, i.e., a hypercube; (2) solution of the encapsulation problem, which is the original governing equations defined in the encapsulation domain; and (3) post-processing of the resulting numerical solution for its statistics, whenever the probability distribution of the epistemic inputs is known, or assumed, a posterior. A notable feature of the approach is that the numerical solution of the encapsulation problem needs to have error control in L1 norm, i.e., point-wise, in the entire encapsulation domain. Though effective and mathematically sound, the requirements of using bounded encapsulation domain and controlling errors in L1 norm can be difficult to achieve in practical simulations. In this paper we extend the capabilities of the methodology proposed in [9]. The new approach achieves two notable extensions. One is in the modeling of the epistemic inputs, where the encapsulation domain is not required to be bounded anymore. Whenever appropriate, unbounded domains can be used to encapsulate the inputs. In practice this could impose less constraints on modeling the inputs, because in many problems there may not be obvious bounds for the epistemic variables. The other extension of the new approach is that the numerical approximation of the encapsulation problem can now be measured in Lp ; p P 1, norm. Naturally this is much easier to achieve in practice than the L1 norm used in the earlier work of [9]. (In fact for semi-bounded and unbounded variables, approximation in L1 norm is not possible.) The ability to use unbounded domain to model the epistemic inputs and to approximate the encapsulation problem in Lp norm makes the new method more flexible in practical computations, as it applies to a much larger set of problems (bounded variables vs. unbounded variables) and allows more flexible numerical treatment (point-wise approximation vs. Lp approximation). In fact the new framework incorporates the earlier one from [9] as a special case, i.e., when the epistemic variables are modeled as bounded variables and a numerical method with point-wise error control is utilized, the new method becomes identical to the earlier one. In term of the analysis of the new method, we present convergence analysis in both weak form and strong form. This is another extension of the earlier work of [9], which only contains convergence result in a weak form. Numerical studies are then presented to demonstrate the properties of the method. 2. Problem setup Let D R‘ ; ‘ ¼ 1; 2; 3, be a physical domain with coordinates x ¼ ðx1 ; . . . ; x‘ Þ and let T > 0 be a real number. We consider the following general stochastic partial differential equation
8 > < v t ðx; t; ZÞ ¼ Lðv Þ; D ð0; T IZ ; Bðv Þ ¼ 0; @D ½0; T IZ ; > : v ¼ v 0; D ft ¼ 0g IZ ;
ð2:1Þ
where L is a (nonlinear) differential operator, B is the boundary condition operator, v 0 is the initial condition, and Z ¼ ðZ 1 ; . . . ; Z d Þ 2 IZ # Rd ; d P 1, are a set of random variables characterizing the random inputs to the governing equation. The solution is therefore a stochastic quantity,
v ðx; t; ZÞ : D ½0; T IZ ! Rnv :
ð2:2Þ
Without loss of generality, hereafter we assume (2.1) is a scalar system with nv ¼ 1. We also make a fundamental assumption that the problem (2.1) is well posed in IZ . Most of the existing studies adopt a probabilistic formulation, which requires the specification of the probability function of the random variables. This translates into the availability of the knowledge of F Z ðsÞ ¼ ProbðZ 6 sÞ, s 2 Rd . (In many practical applications, this is often accomplished by assuming the marginal distributions of each Z i and then assuming mutual independence among all the components.) In this paper, however, we consider the case where the uncertainty is epistemic, in the sense that the distribution function F Z ðsÞ is not completely known. In this case, the standard probabilistic methods do not readily apply.
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
213
Since the focus of our study is on the dependence of the solution on the uncertain inputs Z, hereafter we will suppress the notations of x and t, with the understanding that all statements are meant for any fixed locations in x and t. 3. Methodology We now present our method for solving system (2.1) subject to epistemic uncertain inputs. Similar to the work of [9], the proposed methodology is a three-step procedure that involves identifying the ranges of the uncertain inputs, generating an accurate numerical approximation of the solution to (2.1) within estimated ranges, and post-processing the results. Note that no probability distribution information will be utilized in the solution procedure. 3.1. Range estimation The first task is to identify a range, or a bound, that is sufficiently large such that the ‘‘true’’, and yet unknown, range of the input uncertainty is mostly incorporated. For each variable Z i ; i ¼ 1; . . . ; d, let
IZi ¼ ½ai ; bi ;
ai 6 bi ;
ð3:1Þ
be its (unknown) range. Note that the range can be bounded, semi-bounded, or unbounded. The goal of range estimation is to identify an interval
IX i ¼ ½ai ; bi ;
ð3:2Þ
such that IX i and IZi overlap each other with sufficiently large probability. Note the interval IX i can be semi-bounded or unbounded. This is fundamentally different from the work of [9], where the estimated range IX i is required to be bounded. The estimated range IX i must be sufficiently large such that it encapsulates the ‘‘true’’ range IZi either completely or, if there is any ‘‘truncation’’, the truncation is sufficiently ‘‘small’’. In [9], this requirement is termed the ‘‘overwhelming probability condition’’. In the current setting, the ability to use unbounded interval IX i removes the need for highly accurate estimate of the unknown range IZi , as one can, in principle, always use an unbounded interval to completely encapsulate the unknown range. This could help to alleviate the difficulty in estimating the input ranges, a task that sometimes can be challenging. It should be noted that when estimating the range of each epistemic variables, it is natural to ensure the governing Eq. (2.1) remains well-posed in the estimated range IX i . That is, it certainly makes sense not to let IX i intersect with any parameter values that may violate the physical laws or solvability of the governing equations. For the collection of all the inputs, let IZ be the range of the uncertain epistemic variables Z 2 Rd . Naturally,
IZ # di¼1 IZi :
ð3:3Þ
Note there may be dependence among some of the components of Z, and the real range of Z may be (much) smaller than the tensor products of each range. To encapsulate the range IZ , we define an encapsulation set
IX ¼ di¼1 IX i ¼ di¼1 ½ai ; bi :
ð3:4Þ
which is the Cartesian products of (3.2). Now let
Iþ ¼ IZ [ I X ;
Io ¼ IZ \ IX ;
ð3:5Þ
denoted as super set and common set, respectively, and
I ¼ IZ 4 I X ¼ Iþ n Io ;
ð3:6Þ
be the symmetric difference of IZ and IX , denoted as difference set. The range estimation procedure is deemed satisfactory when
PðZ 2 I Þ 6 d;
ð3:7Þ
for a sufficiently small non-negative real number d P 0. Therefore, the encapsulation set IX encapsulates IZ , the ‘‘true’’ and unknown range of Z, with probability at least 1 d, where d can be made small by enlarging the size of IX . The parameter d can be zero, i.e., IX encapsulates IZ with probability one. One easy way to accomplish this is to ensure IZ # IX , which is relatively easier to enforce in the current framework because IX can be defined as the entire space Rd . 3.2. Encapsulation problem We now define the following encapsulation problem
214
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
8 > < ut ðx; t; XÞ ¼ LðuÞ; D ð0; T IX ; BðuÞ ¼ 0; @D ½0; T IX ; > : u ¼ u0 ; D ft ¼ 0g IX ;
ð3:8Þ
where IX is the bounded hypercube defined in (3.4). This is effectively the same problem (2.1) defined now on the encapsulation set IX that covers the original random parameter set IZ with probability at least 1 d. The new problem is well defined in IX because we have assumed the estimated range of each IX i stays in the range of well-posedness allowed by the governing equation. Since problem (2.1) and (3.8) are exactly the same in the common domain Io , we have the following trivial result,
8s 2 I o :
uð; sÞ ¼ v ð; sÞ;
ð3:9Þ
We remark that for the encapsulation problem (3.8) we do not assign any probability information to the variables X. For the solution of the encapsulation problem (3.8), we focus only on the dependence on the variables X, which now resides in IX # Rd , i.e.,
uðXÞ : IX ! R:
ð3:10Þ
We equip IX with an inner product
ðf ; gÞw ¼
Z
f ðsÞgðsÞwðsÞds;
ð3:11Þ
IX
where wðsÞ P 0 is a weight function, non-vanishing in the interior of IX . We also introduce the weighted Lp norm
kf kLpw ðIX Þ ¼
Z
1=p jf ðsÞjp wðsÞds ;
1 6 p < þ1:
ð3:12Þ
IX
We remark that the choice of the weight function w is entirely a numerical issue. It is used to measure the errors of the numerical approximation and unrelated to the probability information of the uncertain inputs. Let un ðXÞ be a numerical solution to the encapsulation problem (3.8), where n is a discretization parameter that indicates a finer discretization for larger values of n. For example, n can be the highest degree of a polynomial approximation, the total number of discretization elements, etc. A critical requirement for the proposed methodology is the need for the numerical approximation of (3.8) to be accurate in the weighted Lp norm. That is, we require
n , ku un kL
p w ðIX Þ
1:
ð3:13Þ
Ideally the error should converge, i.e., n ! 0 as n ! 1, though convergence is only a mathematical preference, not a practical necessity. This requirement represents another fundamental difference from the method proposed in [9], where the approximate solution un is required to be accurate point-wise in IX , i.e., in L1 norm. The current requirement in the weighted Lp norm is obviously weaker mathematically and easier to accomplish in practice. There exist many numerical approaches that can deliver accurate approximation in the weighted Lp norm. And the choices are problem dependent. Here we will not engage in discussions on the numerical strategies. We assume a satisfactory numerical strategy can always be devised for (3.8) such that the numerical solution is accurate in the weighted Lp norm. We now discuss the way to use un in epistemic uncertainty analysis. 3.3. Epistemic uncertainty analysis The nature of epistemic uncertainty is that it is primarily due to lack of knowledge. And in principle, the true nature of the epistemic uncertainty can be inferred via added information, typically through more measurement and deeper understanding of the system. When un ðXÞ, the polynomial approximation of the true solution uðXÞ, is obtained for (3.8) and accurate in the Lp -norm (3.13), it can serve as an accurate surrogate model. We can then apply various operations on un , instead of u. This is particularly useful when the probability distribution of the variables becomes known a posterior. We can then evaluate the solution statistics of un based on the posterior distribution on Z. Note the operations on un do not require us to solve the governing equations anymore — they can be treated as post-processing steps. And one can repeat this step as much as possible for any kind of posterior distributions. Let qZ ðsÞ ¼ dF Z ðsÞ=ds; s 2 IZ , be the posterior probability distribution of the epistemic uncertain input Z. Then, for example, the mean of the true solution v ðZÞ
l , E½ðv ÞðZÞ ¼
Z
IZ
v ðsÞqZ ðsÞds;
ð3:14Þ
can be approximated by
ln ,
Z Io
un ðsÞqZ ðsÞds:
ð3:15Þ
215
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
Let us define
rðsÞ ¼ IIo ðsÞ
qZ ðsÞ wðsÞ
;
s 2 Iþ ;
ð3:16Þ
where I is the indicator function satisfying, for any set A; IA ðsÞ ¼ 1 for s 2 A and IA ðsÞ ¼ 0 otherwise. Like in [9], where the convergence in mean was established, similar result holds here. Lemma 3.1 (Mean convergence). Assume the solution of (2.1),v ðZÞ, is bounded and let C v ¼ kv kL1 ðIZ Þ . Let un ðXÞ be an approximation to the solution uðXÞ of (3.8) in weighted Lpw norm and denote
n ¼ kun ðXÞ uðXÞkL
p w ðIX Þ
p P 1:
;
ð3:17Þ
Let q be a positive real number satisfying 1=p þ 1=q ¼ 1, and assume the function r defined in (3.16) satisfies, for q < 1 (p > 1),
Z
C r ¼ E½r q1 ¼
Io
rq1 ðsÞqZ ðsÞds < 1:
ð3:18Þ
Then the mean of v in (3.14) and the mean of un in (3.15) satisfy
jl ln j 6 C r1=q n þ C v d;
ð3:19Þ
where for the case of q ¼ 1 (p ¼ 1), we define
C 1=q r
¼ 1.
Proof. We first extend the domains of the definitions for v ; q, and un to Iþ , by following the definitions in (3.5), and define, for s 2 Iþ ,
qþ ðsÞ ¼ IIZ ðsÞqZ ðsÞ;
v þ ðsÞ ¼ II
Z
ðsÞv ðsÞ;
ð3:20Þ
uþn ðsÞ ¼ IIX ðsÞun ðsÞ: Naturally, qþ is a probability density function on Iþ . Then (3.14) can be expressed as
l¼
Z IZ
v ðsÞqZ ðsÞds ¼
Z Iþ
v þ ðsÞqþ ðsÞds;
which can be split into two parts
l¼ ¼
Z Z
Io
Io
Z
v þ ðsÞqþ ðsÞds þ v ðsÞqZ ðsÞds þ
Z
I
I \IZ
v þ ðsÞqþ ðsÞds ¼
Z Io
v ðsÞqZ ðsÞds þ
Z I
v þ ðsÞqþ ðsÞds
v ðsÞqðsÞds:
ð3:21Þ
By using (3.15), we have
l ln ¼
Z Io
ðv ðsÞ un ðsÞÞqZ ðsÞds þ
Z I
v ðsÞqþ ðsÞds ¼
Z Io
ðuðsÞ un ðsÞÞðwðsÞÞ1=p
qZ ðsÞ
ds þ ðwðsÞÞ1=p
Z I \IZ
v ðsÞqZ ðsÞds;
where the property (3.9) has been used. The special case of p ¼ 1 immediately leads to the conclusion. And for p > 1, an exercise of the Hölder inequality leads to
l l 6 n
Z Io
1=p Z juðsÞ un ðsÞjp wðsÞds
Io
qqZ ðsÞ wq=p ðsÞ
1=q ds
þ
Z I \IZ
v ðsÞqZ ðsÞds:
Utilizing (3.7) and (3.18) and the condition of 1=p þ 1=q ¼ 1, the main result (3.19) is established.
h
It is worthwhile to compare this result to that of [9], which states that the error estimate satisfies jl ln j 6 gn þ C v d, where gn is the error of un in L1 -norm. In both results, the second term is the same, indicating the error induced by truncating the true range of the variables. In the current approach, the domain IX can be unbounded, therefore it is easier to avoid this ‘‘truncation error’’ because one does not need to truncate the range. On the other hand, the solution of the encapsulation problem un is required to have a less strict approximation property, in the weighted Lp norm, as opposed to the L1 norm required by [9]. The weight w in the approximation, however, needs to satisfy the condition posed by (3.18), and this manifests itself as the constant in the first term of the error bound (3.19). Additionally, convergence in a strong form can be established. Theorem 3.2 (Strong convergence). Assume the function r, defined in (3.16), satisfies
C r;0 ¼ maxrðsÞ < 1 s
ð3:22Þ
216
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
and the solution of the encapsulation problem un has error in Lpw norm in the form of (3.17). Then 1=p
Eðjv un jp Þ1=p ¼ kv un kLpq ðIþ Þ 6 C r;0 n þ C v d1=p :
ð3:23Þ
Proof
Z jv þ uþn ðsÞjp qþ ðsÞds þ jv þ uþn ðsÞjp qþ ðsÞds Io I Z Z q ðsÞ ðv þ Þp qþ ðsÞds ¼ ju un ðsÞjp wðsÞ Z ds þ ðv þ Þp qþ ðsÞds wðsÞ Io I \IZ Io I \IZ Z p qZ ðsÞ 1=p p p 1=p p þ C ju un ðsÞjp wðsÞds max d 6 C þ C d 6 C þ C d : 6 r;0 n v v v n r;0 s2Io wðsÞ Io
Eðjv un jp Þ ¼
Z
jv þ uþn ðsÞjp qþ ðsÞds ¼ Iþ Z Z ¼ jv un ðsÞjp qþ ðsÞds þ
And this completes the proof.
Z
ð3:24Þ
h
With such a strong convergence in Lpq ðIþ Þ, we immediately obtain the following trivial result. Corollary 3.3 (Weak convergence). Suppose that n and d can be made arbitrarily small. Then, following the same conditions in Theorem 3.2, un converges to v in probability, and consequently also in distribution, with respect to the probability measure qþ defined in (3.20) on the super set Iþ . Naturally, this result implies that the statistical moments of un converges to those of the true solution v, whenever the errors induced by n and d can be arbitrarily small. 3.4. Implementation procedure We now present a quick summary to illustrate the procedure of applying the method to practical problems. Again, we use the setup for the general stochastic system (2.1). Identify the epistemic variables. In (2.1), the variables are denoted as Z ¼ ðZ 1 ; . . . ; Z d Þ; d P 1. (We do not discuss the case of aleatory variables because their treatment is well studied.) For each epistemic variable Z i ; i ¼ 1; . . . ; d, identify a reasonable interval IX i , (3.2), as an estimate of its range. Note this step requires one to utilize any available information, data, or even experience. Define the encapsulation set IX as the tensor products of each interval IX i , as in (3.4). Consequently, this defines the encapsulation problem (3.8). Note that if the intervals IX i are all strictly bounded, then IX is a hypercube and the encapsulation problem becomes identical to that in the earlier work of [9]. Solve the encapsulation problem (3.8), which is the same governing Eq. (2.1) defined in the encapsulated parameter space IX . We now specify a weight function w in IX and seek a numerical solution un that approximates the true solution u in Lpw ; p P 1, norm in IX . This is a standard approximation approach, e.g., finite element methods, and the choices of the weight and the norm vary from problem to problem. Post-process the numerical solution un . Since un is an accurate approximation in a strong norm Lpw , it can be used as a surrogate for the unknown solution u. One can now apply various operations on un to analyze approximately the properties of u. To solve the encapsulation problem (3.8) in the parameter space IX , one can readily borrow several well developed methods such as those based on generalized polynomial chaos, stochastic collocation, etc., all of which offer numerical approximations in a strong norm similar to Lp . Since this is not the focus of the current paper, we refer the interested readers to the book [14] and its large collection of references. 4. Numerical examples In this section, we present numerical examples to support the theoretical analysis. The focus is on examination of the error behavior. We employ two sets of tests. The first one utilizes a scalar equation with a single epistemic variable. Though simple, this examples allows us to thoroughly examine the convergence properties of the numerical implementations, where we employ different variations of encapsulation strategies. The second example is a random diffusion problem with multiple epistemic variables in the diffusivity field. This example allows us to examine the case when the dependence of the variables is unknown. We remark that even though the test problems resemble those in [9], they are now quite different because of the different modeling assumptions on the input epistemic variables. Here we focus on the cases when the epistemic variables are modeled as semi-bounded and unbounded, in order to demonstrate the flexibility of the new framework. Such assumptions are on the epistemic variables are outside the applicability of the framework of [9].
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
217
4.1. Ordinary differential equation Let us consider
dv ¼ Z v ; dt
v ð0Þ ¼ 1;
ð4:1Þ
where the parameter Z > 0 is a random variable representing the input uncertainty. It is an epistemic uncertainty — its probability distribution is unknown a priori. The exact solution is
v ðt; ZÞ ¼ expðZtÞ:
ð4:2Þ
We will use this to examine the error behavior of the numerical methods, by assuming different types of true distributions on Z. The encapsulation problem takes the following form
du ¼ Xu; dt
uð0Þ ¼ 1;
ð4:3Þ
where the range of X; IX , is to be determined based on ones modeling assumptions. Here we consider three difference cases for IX for all possible situations, i.e., IX can be a bounded interval (a,b), a semi-bounded domain ð0; þ1Þ, or an unbounded domain ð1; þ1Þ. Even though this problem admits a simple analytical solution uðt; XÞ ¼ expðXtÞ, we still resort to its numerical solution so that the proposed framework in this paper can be fully examined. Our numerical procedure is a Galerkin method using orthogonal polynomial basis functions. (One is free to use any other suitable numerical methods.) We seek
un ðt; XÞ ¼
n X ^j ðtÞUj ðXÞ; u
ð4:4Þ
j¼0
where fUj ðxÞg are orthogonal polynomials satisfying
Z
Ui ðsÞUj ðsÞwðsÞds ¼ di;j ;
ð4:5Þ
IX
where di;j is the Kronecker delta function. Depending on the range of IX , we utilize different types of orthogonal polynomials. For bounded interval IX ¼ ða; bÞ, we use Legendre polynomials with w const; for semi-bounded interval IX ¼ ð0; þ1Þ, we use Laguerre polynomials with wðsÞ / es ; and for unbounded interval IX ¼ ð1; þ1Þ, we use Hermite polynomials with 2 wðsÞ / es =2 . Upon substituting (4.4) into (4.3) and enforcing the residue to be orthogonal to the linear polynomial space of degree up to n, we obtain the following Galerkin system n X ^j du ^i ; ¼ eij u dt i¼0
^ j ð0Þ ¼ d0;j ; u
ð4:6Þ
R where eij ¼ IX sUi ðsÞUj ðsÞwðsÞds. These are a set of ODEs that can be easily solved. And the resulting numerical solution un converges to u in L2w norm. (See Appendix for proof.) We then assign various distributions to Z a posterior, use these distributions to approximate the statistics of the solution via un directly (without resorting to the numerical simulations), and examine the errors induced by this procedure.
4.1.1. Bounded encapsulation We first assume the exact (and unknown) distribution of Z is an exponential distribution, i.e., qZ ðsÞ expðsÞ. For the encapsulation problem (4.3) we use X 2 ð0; bÞ with various length of b. This is the similar case considered in [9], where the convergence in mean, a weak error measure, was discussed. Here we study the convergence in L2q norm, a strong norm, and examine its behavior based on Theorem 3.2. The results are depicted in Fig. 4.1. It is obvious that the errors converge (exponentially fast) as the order of expansions is increased, before they saturate at certain level. The error saturation level depends on the value of b — a larger value of b induces a smaller error because it causes a smaller ‘‘truncation’’ in the distribution of Z and hence a smaller value d in the second error term in (3.23) of Theorem 3.2. Weak error measure, e.g., error in mean value, behaves similarly. It has been documented in [9] and is not presented here. We further examine the error dependence on b, which is a direct indication of the ‘‘truncation’’ error d. In Fig. 4.2, the strong and weak errors are examined against the parameter b, at a fixed and sufficiently high polynomial order of n ¼ 18. Numerical errors in the Galerkin solver are thus negligible. It can be seen that the errors decrease when b increases, as they should. It is also evident that the strong error converges slower than the weak error with respect to b. This is consistent with the estimate from (3.19) and (3.23), as the strong error scales as d1=p (p ¼ 2 in this case) but weak error scales as d.
218
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
Fig. 4.1. Convergence of the errors in strong L2q norm for linear ODE (4.1) at t ¼ 1, where the encapsulation set is X 2 ð0; bÞ and the exact distribution is exponential qZ ðsÞ expðsÞ.
Fig. 4.2. Convergence of the errors of the strong norm and the weak norm for linear ODE (4.1) at t ¼ 1, where the encapsulation set is X 2 ð0; bÞ and the exact distribution is exponential qZ expðsÞ. The results are for fixed polynomial order of n ¼ 18 and various values of b.
4.1.2. Semi-bounded encapsulation We now employ X 2 ð0; þ1Þ as the encapsulation variable. The corresponds to the practical case where one does not have a reliable estimate of the upper bound of the epistemic variable. We then use Laguerre polynomial with weight wðsÞ / es to approximate the encapsulation problem (4.3). We first assume the true distribution of Z is uniform in (0,1) with qZ ðsÞ ¼ 1. In this case there is no ‘‘truncation’’ error because IX ¼ ð0; þ1Þ completely encapsulates the real variable Z. Error convergence, in both the strong form of L2q and the weak form in mean, are plotted in Fig. 4.3, for different orders of polynomial approximations at a fixed time t ¼ 1. It can be seen that in this case both the strong error and the weak error converge at the same rate – exponential rate – as
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
219
Fig. 4.3. Convergence of errors in the strong form of L2q norm and the weak form of mean, along with the L2w error (n ) of the Galerkin solutions of (4.3), at different polynomial approximation orders n. Here X 2 ð0; 1Þ and the true a posterior distribution of Z is uniform in (0, 1).
the L2w polynomial approximation error of (4.3). This is consistent with the error bounds from Theorems 3.1 and 3.2, as in this case the complete encapsulation of Z by X ensures d ¼ 0. We now assume the true distribution of Z is uniform, but resides in different intervals of (a,b). In Fig. 4.4, the errors are presented from intervals of (0, 1), (5, 6), and (10, 11), in both the L2q form (left) and mean (right). It is obvious that all errors converge (exponentially fast). It also should be remarked that although all intervals have the same size of one, as the interval becomes further away from the origin the errors become larger. This is because in the Laguerre approximation the weight is wðsÞ / es — it achieves the best accuracy close to zero and becomes progressively less accurate away from zero. Consequently, we should emphasize that even though the ability to use semi-bounded and unbounded encapsulation variables in the current framework does alleviate the modeling effort, sometimes significantly, to estimate the ranges of the epistemic variables, it does not imply this can completely eliminate the need to estimate the epistemic variables. When using unbounded domains, the accuracy of the numerical approximation will inevitably be less satisfactory in the region where the weight function approaches zero. It is thus important to have a good, or even just a vague, idea of the true range of the epistemic variables so that one can put the estimated range close to the ‘‘peak’’ of the weight function to ensure more accurate solutions.
Fig. 4.4. Convergence of errors for uniformly distribution Z in different intervals of (a, b), using encapsulation variable X 2 ð0; þ1Þ. Left: L2q error; Right: error in mean.
220
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
4.1.3. Unbounded encapsulation Finally we examine a case where the encapsulation variable is assumed to be X 2 ð1; 1Þ. Correspondingly we employ Hermite polynomial to solve (4.3) numerically. This is a simplified case of the practical situation where one does not have a good understanding of either the upper bound or the lower bound of the epistemic variables. Using an unbounded encapsulation variable would serve as one way to circumvent the difficulty. Both the errors in the strong form (L2q ) and the weak form (mean) are shown in Fig. 4.5, along with the L2w convergence error of the numerical solution of (4.3). Once again we emphasize that in practice it is desirable that one has at least a vague idea of where the true range of the epistemic variable is, and can then ‘‘center’’ the weight function wðsÞ of the Hermite polynomial approximation around that range to achieve more accurate numerical results. 4.2. Homogenous Diffusion Equation We consider the following stochastic diffusion problem in one spatial dimension and multiple random dimensions (d > 1).
@ @x
jðx; ZÞ
@u ðx; ZÞ ¼ f ðx; ZÞ; @x
ðx; ZÞ 2 ð0; 1Þ Rd :
ð4:7Þ
Without loss of generality we impose simple boundary conditions
uð0; ZÞ ¼ 0;
uð1; ZÞ ¼ 0;
Z 2 Rd
and let f ¼ 2 be a constant. For the diffusivity field, we employ the following construction d X 1
jðx; ZÞ ¼ 1 þ r
k¼1
k
2
p2
ð1 þ cosð2pkxÞÞZ k ;
ð4:8Þ
where r > 0 is a parameter to control the variation level of the diffusivity field. Here we choose this fixed form representation to focus on the treatment of the random variables Z k ’s, which are assumed to be epistemic. Note that usually the expansion is obtained by certain decomposition procedure, e.g., Karhunen–Loeve expansion. For non-Gaussian processes, the decomposition usually can not fully determine the probability distributions of Z k ’s. Since the detailed examination of using different encapsulation strategies has been carried out in the previous example, here we mostly focus on dealing with the unknown dependence structure of the epistemic random variables Z k ’s. We set d ¼ 6 and r ¼ 1 in the expansion (4.8). (These choices are rather arbitrary.) We assume that the true (and unknown) distributions of the epistemic variables are: Z 1 betað0; 1; 3; 2Þ, and Z 3 betað1; 0; 1; 1Þ, Z 5 betað0:5; 0:5; 0; 0Þ, and Z 2 ¼ Z 1 Z 5 , Z 4 ¼ ðZ 21 þ 1ÞZ 3 , and Z 6 ¼ Z 5 . Clearly all the variables are bounded, Z 1 2 ð0; 1Þ;
Fig. 4.5. Convergence of the errors of in the strong form of L2q norm and the weak form of mean, along with the L2w error (n ) of the Galerkin solutions of (4.3), at different polynomial approximation orders n. Here the encapsulation set is X 2 ð1; 1Þ and the exact distribution is uniform in (0, 1).
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
221
Fig. 4.6. Convergence of the errors of in both strong form (L2 norm) and in a weak form (in mean) for the diffusion problem. Errors are computed against the high-order (n ¼ 8) Galerkin solution at a fixed spatial location x ¼ 0:5257.
Z 2 2 ð0:5; 0:5Þ, Z 3 2 ð1; 0Þ; Z 4 2 ð2; 0Þ; Z 5 2 ð0:5; 0:5Þ, and Z 6 2 ð0:5; 0:5Þ. Also by construction the variables are dependent. What appears to be a six-dimensional problem is in fact a three-dimensional one, for Z 2 ; Z 4 and Z 6 are functions of the other three variables. In the following epistemic analysis we assume that none of the above information regarding the distributions is unavailable. The only information we have is a conservative estimate of the bounds of the variables. And based on the bound estimation, we employ a straightforward linear transformation and define X 2 IX ¼ ð0; þ1Þ6 to completely encapsulate the epistemic variables Z. Our encapsulation problem is therefore
@ @u ½jðx; XÞ ðx; XÞ ¼ f ðx; XÞ; @x @x
ðx; XÞ 2 ð0; 1Þ ð0; þ1Þ6 ;
ð4:9Þ
where same boundary conditions and source term f hold, and the diffusivity field follows the similar construction of (4.8) with the encapsulation variables X in place of Z. Clearly the encapsulation problem is six-dimensional in the parameter space of X. We employ the six-dimensional Laguerre polynomials, in conjunction with a Galerkin procedure, with weight function Q wðsÞ / 6k¼1 esk to solve the problem. Due to lack of analytical solution, here we employ a sufficiently high-order polynomials solution at n ¼ 8 order as the ‘‘numerical exact’’ solution, against which we compare the solutions obtained by lower-order expansions. Both the strong errors (in L2q ) and weak errors (in mean value) are computed and presented in Fig. 4.6. The errors are computed at a fixed spatial location (x ¼ 0:5257), where the solution is non-trivial. (One can of course take a normed error in the spatial dimension as well.) The usual (exponentially) fast convergence with respect to the approximation order can be clearly observed. No error saturation is present because the encapsulation variables can fully encapsulate the epistemic variables. Various other tests have also be carried out to examine different encapsulation strategies, especially when the encapsulation variables ‘‘truncate’’ the true ranges of the epistemic variables. The results are quite similar to those presented in the previous example, and have also been detailed in the work of [9]. Hence they are not presented here.
5. Summary In this paper we present a computational strategy for epistemic uncertainty analysis. The new method is a notable extension of the work of [9], in the sense that it allows one to use unbounded intervals to encapsulate the true (and unknown) ranges of the epistemic variables. Consequently the new method can employ numerical solutions that are accurate in Lp norms. These feature makes the method more flexible in practice. Theoretical analysis of the method is conducted, where errors in both strong form and weak form are analyzed. Numerical tests are conducted to verify the theory. Though more extensive work is required, and is ongoing, to further examine the methodology, the new method appears to be a fairly effective tool for epistemic uncertainty analysis.
222
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
Appendix: Convergence of Galerkin solution for ODE We now present the convergence proof for the orthogonal polynomial Galerkin solution for the ordinary differential Eq. (4.3). In fact, we consider a slightly more general form of the problem
du ¼ aðXÞu; dt
0 < t 6 T;
uðt ¼ 0Þ ¼ 1;
ð5:1Þ
where the coefficient a is bounded away from infinity
1 < amin 6 a 6 amax < þ1:
ð5:2Þ
Let IX be the range of X and fUj g be a set of orthogonal polynomials defined on IX
Z
Uj ðsÞUk ðsÞwðsÞds ¼ djk ;
ð5:3Þ
IX
R where w > 0 is a weight function inducing an inner product ðf ; gÞw ¼ f ðsÞgðsÞwðsÞds and norm kf kw ¼ ðf ; f Þ1=2 w . Here and m after, the space Hw ðIX Þ denotes the usual weighted Sobolev spaces equipped with the norm kf kHmw . In practice, it is often convenient to map IX to a standard interval such as ð1; þ1Þ, ð0; þ1Þ, or ð1; 1Þ, where standard forms of orthogonal polynomial can be readily applied. Hereafter we will not explicitly specify the interval. In Gelerkin method we seek an expansion
v N ðt; XÞ ¼
N X
v^ j ðtÞUj ðXÞ
ð5:4Þ
j¼0
such that (5.1) is satisfied in a weak form, i.e., the residue is orthogonal to the polynomial space where the solution is constructed. This leads to, for all k ¼ 0; . . . ; N,
dv N ¼ 0: þ av N ; Uk dt w Straightforward derivation then gives us, for all k ¼ 0; . . . ; N, N X dv^ k ¼ ajk v^ j ; dt j¼0
ð5:5Þ
where
ajk ¼
Z
aðsÞUj ðsÞUk ðsÞwðsÞds
ð5:6Þ
^ k ðt ¼ 0Þ ¼ dk0 . and the initial conditions are v In order to prove the convergence of the Galerkin solution v N to the solution u of (5.1), it is convenient to consider an intermediate projection. More precisely, for uðt; XÞ 2 L2w ðIX Þ for any t, its orthogonal projection P N u exists and is defined as
PN uðt; XÞ ,
N X
^ j ðtÞUj ðXÞ; u
^ j ¼ ðu; Uj Þw : u
ð5:7Þ
j¼0
Obviously, P N u is optimal in the k kw norm and limN!1 P N u ¼ u. Assume that for some m > 1=2,
kðPN u uÞðtÞkw 6 CN m kuðtÞkHmw ðIX Þ
8 t > 0;
ð5:8Þ
where the constant C does not depend on N. Note that the rate of convergence depends on the regularity of u and the type of orthogonal polynomials fUj g. Theorem 5.1 (Convergence). Suppose that the assumption (5.8) holds. The Galerkin solution v N from (5.4)–(5.6) converges to the solution u of the Eq. (5.1). Moreover, if the solution uðt; XÞ belongs to the weighted Sobolev space Hm w ðI X Þ for any t, there exists a constant C, independent of N, such that
maxkðv N uÞðtÞkw 6 CNmþ1=2 maxkuðtÞkHmw ðIX Þ :
06t6T
06t6T
ð5:9Þ
^ j g given by (5.7) satisfy the following equations, Proof. It is easy to verify that the coefficients fu 1 X ^k du ^j ; ¼ ajk u dt j¼0
k ¼ 0; 1; . . . ; 1:
ð5:10Þ
223
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
Let us define
^k ; ^ek ¼ v^ k u
k ¼ 0; . . . ; N:
ð5:11Þ
We then have N 1 X X d^ek ^j : ajk u ¼ ajk ^ej þ dt j¼0 j>N
Upon multiplying both sides by ^ ek and taking a summation over k, we obtain N N X N N X 1 X X X d^e ^ek k ¼ ^ek ajk ^ej þ ^ek ajk u ^j : dt k¼0 k¼0 j¼0 k¼0 j>N
Let eN ðt; XÞ ¼ v N P N u. Naturally we have keN kw ¼
e ¼ ð^e0 ; . . . ; ^eN ÞT ;
PN
^2 k¼0 ek .
We now adopt vector–matrix notation by defining
A ¼ ðajk Þ06j;k6N :
The above equation can then be written as
1 d kek2 ¼ eT Ae þ eT r; 2 dt
ð5:12Þ
where the norm is the vector 2-norm and kek ¼ keN kw and r ¼ ð^r0 ; . . . ; ^r N ÞT with the entries
^rk ¼
1 X ^j ; ajk u
k ¼ 0; . . . ; N:
ð5:13Þ
j>N
Lemma 5.2. Following the definitions of e and A and the condition (5.2)
eT Ae P amin kek2 :
ð5:14Þ
Proof. Consider an arbitrary vector e and construct a function
qðsÞ ¼ ^e0 þ ^e1 U1 ðsÞ þ þ ^eN UN ðsÞ: Obviously qðsÞ is a polynomial and satisfies kqkw ¼ kek.
eT Ae ¼
Z N X N N X N Z X X ^ej ajk ^ek ¼ aðsÞ^ej ^ek Uj ðsÞUk ðsÞwðsÞds ¼ aðsÞq2 ðsÞwðsÞds P amin kqk2w : j¼0 k¼0
ð5:15Þ
j¼0 k¼0
h Substituting this result into (5.12), we obtain
d kek 6 amin kek þ krk: dt
ð5:16Þ
We now analyze the term krk. Using the definition (5.13), we derive
^rk ¼
Z Z 1 1 Z 1 X X X ^j ¼ ajk u aðsÞUj ðsÞUk ðsÞwðsÞdsu^ j ¼ aðsÞ u^ j Uj ðsÞUk ðsÞwðsÞds ¼ aðsÞðu PN uÞUk ðsÞwðsÞds: j>N
j>N
j>N
Then, it follows that
krk ¼
N X ^r2k k¼0
!1=2 ¼
N Z X k¼0
2 !1=2 6 ðu PN uÞaðsÞUk ðsÞwðsÞds
pffiffiffiffiffiffiffiffiffiffiffiffiffi 6 amax N þ 1ku PN ukw ; , bðN; tÞ;
N X ku PN uk2w kaUk k2w
!1=2
k¼0
ð5:17Þ
where the Cauchy–Schwarz inequality and the orthonormal property of the basis are used. Note that for m > 1=2; bðN; tÞ ! 0 as N ! 1 for all t due to (5.8). And from (5.16 and 5.17), we obtain
d kek 6 amin kek þ bðN; tÞ: dt
ð5:18Þ
224
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
The Gronwall inequality then leads to
keðtÞk 6 eamin t keð0Þk þ
Z
t
eamin ðtsÞ bðN; sÞds:
ð5:19Þ
0
An application of the triangle inequality implies
kv N ukw 6 kv N PN ukw þ kPN u ukw and the convergence theorem is then established for N ! 1 by using (5.8), (5.17) and keð0Þk ¼ 0.
h
It is trivial to see that the same convergence result holds true for the case of keð0Þk – 0 but converges at the same rate of the projection error ku PN ukw . This corresponds to the case of (5.1) with a random initial condition u0 ðXÞ that is approximated by its orthogonal projection PN u0 . Correspondingly the initial conditions of the system (5.5) becomes v^ k ðt ¼ 0Þ ¼ ðu0 ; Uk Þw . This is the typical procedure in practice. References [1] I. Babuska, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal. 45 (3) (2007) 1005–1034. [2] I. Babuska, R. Tempone, G.E. Zouraris, Galerkin finite element approximations of stochastic elliptic differential equations, SIAM J. Numer. Anal. 42 (2004) 800–825. [3] D. Dubois, H. Prade, Possibility Theory: An Approach to Computerized Processing of Uncertainty, Plenum, New York, 1998. [4] M. Eldred, Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design, in: 50th AIAA/ ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, vol. AIAA-2009-2249, 2009. [5] M. Eldred, L. Swilder, G. Tang, Mixed aleatory-epistemic uncertainty quantification with stochastic expansions and optimization-based interval estimation, Reliab. Eng. Syst. Saf. 96 (9) (2011) 1092–1113. [6] R.G. Ghanem, P. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer-Verlag, 1991. [7] J.C. Helton, J.D. Johnson, W.L. Oberkampf, C.J. Sallaberry, Representation of analysis results involving aleatory and epistemic uncertainty, Technical Report 4379, Sandia National Laboratories, 2008. [8] J.C. Helton, J.D. Johnson, W.L. Oberkampf, C.B. Storlie, A sampling-based computational strategy for the representation of epistemic uncertainty in model predictions with evidence theory, Technical Report 5557, Sandia National Laboratories, 2006. [9] J. Jakeman, M. Eldred, D. Xiu, Numerical approach for quantification of epistemic uncertainty, J. Comput. Phys. 229 (2010) 4648–4663. [10] L. Swiler, T. Paez, R. Mayes, M. Eldred. Epistemic uncertainty in the calculation of margins, in: AIAA Structures, Structural Dynamics, and Materials Conference, Palm Springs CA, May 2009. [11] L.P. Swiler, T.L. Paez, R.L. Mayes, Epistemic uncertainty quantification tutorial, in: IMAC XXVII Conference and Exposition on Structural Dynamics, Orlando, FL, Feburary 2009. [12] D. Xiu, Efficient collocational approach for parametric uncertainty analysis, Commun. Comput. Phys. 2 (2) (2007) 293–309. [13] D. Xiu, Fast numerical methods for stochastic computations: a review, Commun. Comput. Phys. 5 (2009) 242–272. [14] D. Xiu, Numerical Methods for Stochastic Computations, Princeton Univeristy Press, Princeton, New Jersey, 2010. [15] D. Xiu, J.S. Hesthaven, High-order collocation methods for differential equations with random inputs, SIAM J. Sci. Comput. 27 (3) (2005) 1118–1139. [16] D. Xiu, G.E. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput. 24 (2) (2002) 619–644.
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
A flexible numerical approach for quantification of epistemic uncertainty Xiaoxiao Chen a, Eun-Jae Park b,1, Dongbin Xiu a,c,⇑,2 a
Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA Department of Computational Science & Engineering, Yonsei University, Seoul 120-749, Republic of Korea c Department of Mathematics, and Scientific Computing and Imagining Institute University of Utah, Salt Lake City, UT 84112, USA b
a r t i c l e
i n f o
Article history: Received 2 December 2011 Received in revised form 15 January 2013 Accepted 19 January 2013 Available online 4 February 2013 Keywords: Uncertainty quantification Epistemic uncertainty Generalized polynomial chaos Encapsulation problem
a b s t r a c t In the field of uncertainty quantification (UQ), epistemic uncertainty often refers to the kind of uncertainty whose complete probabilistic description is not available, largely due to our lack of knowledge about the uncertainty. Quantification of the impacts of epistemic uncertainty is naturally difficult, because most of the existing stochastic tools rely on the specification of the probability distributions and thus do not readily apply to epistemic uncertainty. And there have been few studies and methods to deal with epistemic uncertainty. A recent work can be found in [J. Jakeman, M. Eldred, D. Xiu, Numerical approach for quantification of epistemic uncertainty, J. Comput. Phys. 229 (2010) 4648–4663], where a framework for numerical treatment of epistemic uncertainty was proposed. The method is based on solving an encapsulation problem, without using any probability information, in a hypercube that encapsulates the unknown epistemic probability space. If more probabilistic information about the epistemic variables is known a posteriori, the solution statistics can then be evaluated at post-process steps. In this paper, we present a new method, similar to that of Jakeman et al. but significantly extending its capabilities. Most notably, the new method (1) does not require the encapsulation problem to be in a bounded domain such as a hypercube; (2) does not require the solution of the encapsulation problem to converge point-wise. In the current formulation, the encapsulation problem could reside in an unbounded domain, and more importantly, its numerical approximation could be sought in Lp norm. These features thus make the new approach more flexible and amicable to practical implementation. Both the mathematical framework and numerical analysis are presented to demonstrate the effectiveness of the new approach. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Extensive research efforts have been devoted to uncertainty quantification (UQ), especially for large and complex systems. One of the key issues in UQ is to design efficient numerical algorithms to propagate uncertainty from system inputs to outputs. To this end, the most widely used approach adopts probabilistic framework, where the uncertainty inputs are modeled as random variables/processes. One condition naturally required by most probabilistic methods is the knowledge ⇑ Corresponding author. Tel.: +1 765 496 2846. E-mail addresses: [email protected] (E.-J. Park), [email protected]; [email protected] (D. Xiu). The work of this author was supported in part by the WCU program through the Korea Science and Engineering Foundation funded by the Ministry of Education, Science, and Technology R31-2008-000-10049-0. 2 The work of D. Xiu was supported by AFOSR, DOE, NNSA, and NSF; also by NRF R31-2008-000-10049-0 through Department of Computational Science & Engineering, Yonsei University, Seoul 120-749, Korea. Tel.: +1 202 642 4808. 1
0021-9991/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.01.018
212
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
of the probability distribution of the uncertain inputs. In the field of UQ, this kind of uncertainty is sometimes referred to as aleatory uncertainty — the kind of uncertainty whose complete probabilistic specification is available. In practice, however, it is often very difficult, if not impossible, to have sufficient information to fully prescribe the probability distribution of all the random inputs. And in such cases, one is facing epistemic uncertainty, the kind of uncertainty whose probability distribution is not known completely. Until recently, most uncertainty analysis has focused on aleatory uncertainty. Numerous stochastic methods have been developed to provide accurate and efficient simulations for this form of uncertainty. For example, the stochastic methods based on generalized polynomial chaos (gPC) [16], an extension of the traditional polynomial chaos [6], have demonstrated their efficiency in practice and are widely used for many problems. Their implementations typically follow either stochastic Galerkin (SG) [2,6,16] or stochastic collocation (SC) [1,4,12,15], often termed as intrusive methods or non-intrusive methods, respectively. For a detailed review on the methods, see [13,14]. Numerical study of epistemic uncertainty is more difficult because of the lack of sufficient probabilistic information. And one naturally can not readily adopt probabilistic approaches. Some of the existing approaches include evidence theory [8], possibility theory [3] and interval analysis [7,10]. These method have their own advantages, though most do not address efficient numerical implementations. More recent studies can be found in [5,11]. This paper is based on a more recent work of [9], where a general numerical approach was proposed for epistemic uncertainty analysis. The method consists of three components: (1) encapsulation of the epistemic uncertain inputs by a bounded domain, i.e., a hypercube; (2) solution of the encapsulation problem, which is the original governing equations defined in the encapsulation domain; and (3) post-processing of the resulting numerical solution for its statistics, whenever the probability distribution of the epistemic inputs is known, or assumed, a posterior. A notable feature of the approach is that the numerical solution of the encapsulation problem needs to have error control in L1 norm, i.e., point-wise, in the entire encapsulation domain. Though effective and mathematically sound, the requirements of using bounded encapsulation domain and controlling errors in L1 norm can be difficult to achieve in practical simulations. In this paper we extend the capabilities of the methodology proposed in [9]. The new approach achieves two notable extensions. One is in the modeling of the epistemic inputs, where the encapsulation domain is not required to be bounded anymore. Whenever appropriate, unbounded domains can be used to encapsulate the inputs. In practice this could impose less constraints on modeling the inputs, because in many problems there may not be obvious bounds for the epistemic variables. The other extension of the new approach is that the numerical approximation of the encapsulation problem can now be measured in Lp ; p P 1, norm. Naturally this is much easier to achieve in practice than the L1 norm used in the earlier work of [9]. (In fact for semi-bounded and unbounded variables, approximation in L1 norm is not possible.) The ability to use unbounded domain to model the epistemic inputs and to approximate the encapsulation problem in Lp norm makes the new method more flexible in practical computations, as it applies to a much larger set of problems (bounded variables vs. unbounded variables) and allows more flexible numerical treatment (point-wise approximation vs. Lp approximation). In fact the new framework incorporates the earlier one from [9] as a special case, i.e., when the epistemic variables are modeled as bounded variables and a numerical method with point-wise error control is utilized, the new method becomes identical to the earlier one. In term of the analysis of the new method, we present convergence analysis in both weak form and strong form. This is another extension of the earlier work of [9], which only contains convergence result in a weak form. Numerical studies are then presented to demonstrate the properties of the method. 2. Problem setup Let D R‘ ; ‘ ¼ 1; 2; 3, be a physical domain with coordinates x ¼ ðx1 ; . . . ; x‘ Þ and let T > 0 be a real number. We consider the following general stochastic partial differential equation
8 > < v t ðx; t; ZÞ ¼ Lðv Þ; D ð0; T IZ ; Bðv Þ ¼ 0; @D ½0; T IZ ; > : v ¼ v 0; D ft ¼ 0g IZ ;
ð2:1Þ
where L is a (nonlinear) differential operator, B is the boundary condition operator, v 0 is the initial condition, and Z ¼ ðZ 1 ; . . . ; Z d Þ 2 IZ # Rd ; d P 1, are a set of random variables characterizing the random inputs to the governing equation. The solution is therefore a stochastic quantity,
v ðx; t; ZÞ : D ½0; T IZ ! Rnv :
ð2:2Þ
Without loss of generality, hereafter we assume (2.1) is a scalar system with nv ¼ 1. We also make a fundamental assumption that the problem (2.1) is well posed in IZ . Most of the existing studies adopt a probabilistic formulation, which requires the specification of the probability function of the random variables. This translates into the availability of the knowledge of F Z ðsÞ ¼ ProbðZ 6 sÞ, s 2 Rd . (In many practical applications, this is often accomplished by assuming the marginal distributions of each Z i and then assuming mutual independence among all the components.) In this paper, however, we consider the case where the uncertainty is epistemic, in the sense that the distribution function F Z ðsÞ is not completely known. In this case, the standard probabilistic methods do not readily apply.
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
213
Since the focus of our study is on the dependence of the solution on the uncertain inputs Z, hereafter we will suppress the notations of x and t, with the understanding that all statements are meant for any fixed locations in x and t. 3. Methodology We now present our method for solving system (2.1) subject to epistemic uncertain inputs. Similar to the work of [9], the proposed methodology is a three-step procedure that involves identifying the ranges of the uncertain inputs, generating an accurate numerical approximation of the solution to (2.1) within estimated ranges, and post-processing the results. Note that no probability distribution information will be utilized in the solution procedure. 3.1. Range estimation The first task is to identify a range, or a bound, that is sufficiently large such that the ‘‘true’’, and yet unknown, range of the input uncertainty is mostly incorporated. For each variable Z i ; i ¼ 1; . . . ; d, let
IZi ¼ ½ai ; bi ;
ai 6 bi ;
ð3:1Þ
be its (unknown) range. Note that the range can be bounded, semi-bounded, or unbounded. The goal of range estimation is to identify an interval
IX i ¼ ½ai ; bi ;
ð3:2Þ
such that IX i and IZi overlap each other with sufficiently large probability. Note the interval IX i can be semi-bounded or unbounded. This is fundamentally different from the work of [9], where the estimated range IX i is required to be bounded. The estimated range IX i must be sufficiently large such that it encapsulates the ‘‘true’’ range IZi either completely or, if there is any ‘‘truncation’’, the truncation is sufficiently ‘‘small’’. In [9], this requirement is termed the ‘‘overwhelming probability condition’’. In the current setting, the ability to use unbounded interval IX i removes the need for highly accurate estimate of the unknown range IZi , as one can, in principle, always use an unbounded interval to completely encapsulate the unknown range. This could help to alleviate the difficulty in estimating the input ranges, a task that sometimes can be challenging. It should be noted that when estimating the range of each epistemic variables, it is natural to ensure the governing Eq. (2.1) remains well-posed in the estimated range IX i . That is, it certainly makes sense not to let IX i intersect with any parameter values that may violate the physical laws or solvability of the governing equations. For the collection of all the inputs, let IZ be the range of the uncertain epistemic variables Z 2 Rd . Naturally,
IZ # di¼1 IZi :
ð3:3Þ
Note there may be dependence among some of the components of Z, and the real range of Z may be (much) smaller than the tensor products of each range. To encapsulate the range IZ , we define an encapsulation set
IX ¼ di¼1 IX i ¼ di¼1 ½ai ; bi :
ð3:4Þ
which is the Cartesian products of (3.2). Now let
Iþ ¼ IZ [ I X ;
Io ¼ IZ \ IX ;
ð3:5Þ
denoted as super set and common set, respectively, and
I ¼ IZ 4 I X ¼ Iþ n Io ;
ð3:6Þ
be the symmetric difference of IZ and IX , denoted as difference set. The range estimation procedure is deemed satisfactory when
PðZ 2 I Þ 6 d;
ð3:7Þ
for a sufficiently small non-negative real number d P 0. Therefore, the encapsulation set IX encapsulates IZ , the ‘‘true’’ and unknown range of Z, with probability at least 1 d, where d can be made small by enlarging the size of IX . The parameter d can be zero, i.e., IX encapsulates IZ with probability one. One easy way to accomplish this is to ensure IZ # IX , which is relatively easier to enforce in the current framework because IX can be defined as the entire space Rd . 3.2. Encapsulation problem We now define the following encapsulation problem
214
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
8 > < ut ðx; t; XÞ ¼ LðuÞ; D ð0; T IX ; BðuÞ ¼ 0; @D ½0; T IX ; > : u ¼ u0 ; D ft ¼ 0g IX ;
ð3:8Þ
where IX is the bounded hypercube defined in (3.4). This is effectively the same problem (2.1) defined now on the encapsulation set IX that covers the original random parameter set IZ with probability at least 1 d. The new problem is well defined in IX because we have assumed the estimated range of each IX i stays in the range of well-posedness allowed by the governing equation. Since problem (2.1) and (3.8) are exactly the same in the common domain Io , we have the following trivial result,
8s 2 I o :
uð; sÞ ¼ v ð; sÞ;
ð3:9Þ
We remark that for the encapsulation problem (3.8) we do not assign any probability information to the variables X. For the solution of the encapsulation problem (3.8), we focus only on the dependence on the variables X, which now resides in IX # Rd , i.e.,
uðXÞ : IX ! R:
ð3:10Þ
We equip IX with an inner product
ðf ; gÞw ¼
Z
f ðsÞgðsÞwðsÞds;
ð3:11Þ
IX
where wðsÞ P 0 is a weight function, non-vanishing in the interior of IX . We also introduce the weighted Lp norm
kf kLpw ðIX Þ ¼
Z
1=p jf ðsÞjp wðsÞds ;
1 6 p < þ1:
ð3:12Þ
IX
We remark that the choice of the weight function w is entirely a numerical issue. It is used to measure the errors of the numerical approximation and unrelated to the probability information of the uncertain inputs. Let un ðXÞ be a numerical solution to the encapsulation problem (3.8), where n is a discretization parameter that indicates a finer discretization for larger values of n. For example, n can be the highest degree of a polynomial approximation, the total number of discretization elements, etc. A critical requirement for the proposed methodology is the need for the numerical approximation of (3.8) to be accurate in the weighted Lp norm. That is, we require
n , ku un kL
p w ðIX Þ
1:
ð3:13Þ
Ideally the error should converge, i.e., n ! 0 as n ! 1, though convergence is only a mathematical preference, not a practical necessity. This requirement represents another fundamental difference from the method proposed in [9], where the approximate solution un is required to be accurate point-wise in IX , i.e., in L1 norm. The current requirement in the weighted Lp norm is obviously weaker mathematically and easier to accomplish in practice. There exist many numerical approaches that can deliver accurate approximation in the weighted Lp norm. And the choices are problem dependent. Here we will not engage in discussions on the numerical strategies. We assume a satisfactory numerical strategy can always be devised for (3.8) such that the numerical solution is accurate in the weighted Lp norm. We now discuss the way to use un in epistemic uncertainty analysis. 3.3. Epistemic uncertainty analysis The nature of epistemic uncertainty is that it is primarily due to lack of knowledge. And in principle, the true nature of the epistemic uncertainty can be inferred via added information, typically through more measurement and deeper understanding of the system. When un ðXÞ, the polynomial approximation of the true solution uðXÞ, is obtained for (3.8) and accurate in the Lp -norm (3.13), it can serve as an accurate surrogate model. We can then apply various operations on un , instead of u. This is particularly useful when the probability distribution of the variables becomes known a posterior. We can then evaluate the solution statistics of un based on the posterior distribution on Z. Note the operations on un do not require us to solve the governing equations anymore — they can be treated as post-processing steps. And one can repeat this step as much as possible for any kind of posterior distributions. Let qZ ðsÞ ¼ dF Z ðsÞ=ds; s 2 IZ , be the posterior probability distribution of the epistemic uncertain input Z. Then, for example, the mean of the true solution v ðZÞ
l , E½ðv ÞðZÞ ¼
Z
IZ
v ðsÞqZ ðsÞds;
ð3:14Þ
can be approximated by
ln ,
Z Io
un ðsÞqZ ðsÞds:
ð3:15Þ
215
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
Let us define
rðsÞ ¼ IIo ðsÞ
qZ ðsÞ wðsÞ
;
s 2 Iþ ;
ð3:16Þ
where I is the indicator function satisfying, for any set A; IA ðsÞ ¼ 1 for s 2 A and IA ðsÞ ¼ 0 otherwise. Like in [9], where the convergence in mean was established, similar result holds here. Lemma 3.1 (Mean convergence). Assume the solution of (2.1),v ðZÞ, is bounded and let C v ¼ kv kL1 ðIZ Þ . Let un ðXÞ be an approximation to the solution uðXÞ of (3.8) in weighted Lpw norm and denote
n ¼ kun ðXÞ uðXÞkL
p w ðIX Þ
p P 1:
;
ð3:17Þ
Let q be a positive real number satisfying 1=p þ 1=q ¼ 1, and assume the function r defined in (3.16) satisfies, for q < 1 (p > 1),
Z
C r ¼ E½r q1 ¼
Io
rq1 ðsÞqZ ðsÞds < 1:
ð3:18Þ
Then the mean of v in (3.14) and the mean of un in (3.15) satisfy
jl ln j 6 C r1=q n þ C v d;
ð3:19Þ
where for the case of q ¼ 1 (p ¼ 1), we define
C 1=q r
¼ 1.
Proof. We first extend the domains of the definitions for v ; q, and un to Iþ , by following the definitions in (3.5), and define, for s 2 Iþ ,
qþ ðsÞ ¼ IIZ ðsÞqZ ðsÞ;
v þ ðsÞ ¼ II
Z
ðsÞv ðsÞ;
ð3:20Þ
uþn ðsÞ ¼ IIX ðsÞun ðsÞ: Naturally, qþ is a probability density function on Iþ . Then (3.14) can be expressed as
l¼
Z IZ
v ðsÞqZ ðsÞds ¼
Z Iþ
v þ ðsÞqþ ðsÞds;
which can be split into two parts
l¼ ¼
Z Z
Io
Io
Z
v þ ðsÞqþ ðsÞds þ v ðsÞqZ ðsÞds þ
Z
I
I \IZ
v þ ðsÞqþ ðsÞds ¼
Z Io
v ðsÞqZ ðsÞds þ
Z I
v þ ðsÞqþ ðsÞds
v ðsÞqðsÞds:
ð3:21Þ
By using (3.15), we have
l ln ¼
Z Io
ðv ðsÞ un ðsÞÞqZ ðsÞds þ
Z I
v ðsÞqþ ðsÞds ¼
Z Io
ðuðsÞ un ðsÞÞðwðsÞÞ1=p
qZ ðsÞ
ds þ ðwðsÞÞ1=p
Z I \IZ
v ðsÞqZ ðsÞds;
where the property (3.9) has been used. The special case of p ¼ 1 immediately leads to the conclusion. And for p > 1, an exercise of the Hölder inequality leads to
l l 6 n
Z Io
1=p Z juðsÞ un ðsÞjp wðsÞds
Io
qqZ ðsÞ wq=p ðsÞ
1=q ds
þ
Z I \IZ
v ðsÞqZ ðsÞds:
Utilizing (3.7) and (3.18) and the condition of 1=p þ 1=q ¼ 1, the main result (3.19) is established.
h
It is worthwhile to compare this result to that of [9], which states that the error estimate satisfies jl ln j 6 gn þ C v d, where gn is the error of un in L1 -norm. In both results, the second term is the same, indicating the error induced by truncating the true range of the variables. In the current approach, the domain IX can be unbounded, therefore it is easier to avoid this ‘‘truncation error’’ because one does not need to truncate the range. On the other hand, the solution of the encapsulation problem un is required to have a less strict approximation property, in the weighted Lp norm, as opposed to the L1 norm required by [9]. The weight w in the approximation, however, needs to satisfy the condition posed by (3.18), and this manifests itself as the constant in the first term of the error bound (3.19). Additionally, convergence in a strong form can be established. Theorem 3.2 (Strong convergence). Assume the function r, defined in (3.16), satisfies
C r;0 ¼ maxrðsÞ < 1 s
ð3:22Þ
216
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
and the solution of the encapsulation problem un has error in Lpw norm in the form of (3.17). Then 1=p
Eðjv un jp Þ1=p ¼ kv un kLpq ðIþ Þ 6 C r;0 n þ C v d1=p :
ð3:23Þ
Proof
Z jv þ uþn ðsÞjp qþ ðsÞds þ jv þ uþn ðsÞjp qþ ðsÞds Io I Z Z q ðsÞ ðv þ Þp qþ ðsÞds ¼ ju un ðsÞjp wðsÞ Z ds þ ðv þ Þp qþ ðsÞds wðsÞ Io I \IZ Io I \IZ Z p qZ ðsÞ 1=p p p 1=p p þ C ju un ðsÞjp wðsÞds max d 6 C þ C d 6 C þ C d : 6 r;0 n v v v n r;0 s2Io wðsÞ Io
Eðjv un jp Þ ¼
Z
jv þ uþn ðsÞjp qþ ðsÞds ¼ Iþ Z Z ¼ jv un ðsÞjp qþ ðsÞds þ
And this completes the proof.
Z
ð3:24Þ
h
With such a strong convergence in Lpq ðIþ Þ, we immediately obtain the following trivial result. Corollary 3.3 (Weak convergence). Suppose that n and d can be made arbitrarily small. Then, following the same conditions in Theorem 3.2, un converges to v in probability, and consequently also in distribution, with respect to the probability measure qþ defined in (3.20) on the super set Iþ . Naturally, this result implies that the statistical moments of un converges to those of the true solution v, whenever the errors induced by n and d can be arbitrarily small. 3.4. Implementation procedure We now present a quick summary to illustrate the procedure of applying the method to practical problems. Again, we use the setup for the general stochastic system (2.1). Identify the epistemic variables. In (2.1), the variables are denoted as Z ¼ ðZ 1 ; . . . ; Z d Þ; d P 1. (We do not discuss the case of aleatory variables because their treatment is well studied.) For each epistemic variable Z i ; i ¼ 1; . . . ; d, identify a reasonable interval IX i , (3.2), as an estimate of its range. Note this step requires one to utilize any available information, data, or even experience. Define the encapsulation set IX as the tensor products of each interval IX i , as in (3.4). Consequently, this defines the encapsulation problem (3.8). Note that if the intervals IX i are all strictly bounded, then IX is a hypercube and the encapsulation problem becomes identical to that in the earlier work of [9]. Solve the encapsulation problem (3.8), which is the same governing Eq. (2.1) defined in the encapsulated parameter space IX . We now specify a weight function w in IX and seek a numerical solution un that approximates the true solution u in Lpw ; p P 1, norm in IX . This is a standard approximation approach, e.g., finite element methods, and the choices of the weight and the norm vary from problem to problem. Post-process the numerical solution un . Since un is an accurate approximation in a strong norm Lpw , it can be used as a surrogate for the unknown solution u. One can now apply various operations on un to analyze approximately the properties of u. To solve the encapsulation problem (3.8) in the parameter space IX , one can readily borrow several well developed methods such as those based on generalized polynomial chaos, stochastic collocation, etc., all of which offer numerical approximations in a strong norm similar to Lp . Since this is not the focus of the current paper, we refer the interested readers to the book [14] and its large collection of references. 4. Numerical examples In this section, we present numerical examples to support the theoretical analysis. The focus is on examination of the error behavior. We employ two sets of tests. The first one utilizes a scalar equation with a single epistemic variable. Though simple, this examples allows us to thoroughly examine the convergence properties of the numerical implementations, where we employ different variations of encapsulation strategies. The second example is a random diffusion problem with multiple epistemic variables in the diffusivity field. This example allows us to examine the case when the dependence of the variables is unknown. We remark that even though the test problems resemble those in [9], they are now quite different because of the different modeling assumptions on the input epistemic variables. Here we focus on the cases when the epistemic variables are modeled as semi-bounded and unbounded, in order to demonstrate the flexibility of the new framework. Such assumptions are on the epistemic variables are outside the applicability of the framework of [9].
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
217
4.1. Ordinary differential equation Let us consider
dv ¼ Z v ; dt
v ð0Þ ¼ 1;
ð4:1Þ
where the parameter Z > 0 is a random variable representing the input uncertainty. It is an epistemic uncertainty — its probability distribution is unknown a priori. The exact solution is
v ðt; ZÞ ¼ expðZtÞ:
ð4:2Þ
We will use this to examine the error behavior of the numerical methods, by assuming different types of true distributions on Z. The encapsulation problem takes the following form
du ¼ Xu; dt
uð0Þ ¼ 1;
ð4:3Þ
where the range of X; IX , is to be determined based on ones modeling assumptions. Here we consider three difference cases for IX for all possible situations, i.e., IX can be a bounded interval (a,b), a semi-bounded domain ð0; þ1Þ, or an unbounded domain ð1; þ1Þ. Even though this problem admits a simple analytical solution uðt; XÞ ¼ expðXtÞ, we still resort to its numerical solution so that the proposed framework in this paper can be fully examined. Our numerical procedure is a Galerkin method using orthogonal polynomial basis functions. (One is free to use any other suitable numerical methods.) We seek
un ðt; XÞ ¼
n X ^j ðtÞUj ðXÞ; u
ð4:4Þ
j¼0
where fUj ðxÞg are orthogonal polynomials satisfying
Z
Ui ðsÞUj ðsÞwðsÞds ¼ di;j ;
ð4:5Þ
IX
where di;j is the Kronecker delta function. Depending on the range of IX , we utilize different types of orthogonal polynomials. For bounded interval IX ¼ ða; bÞ, we use Legendre polynomials with w const; for semi-bounded interval IX ¼ ð0; þ1Þ, we use Laguerre polynomials with wðsÞ / es ; and for unbounded interval IX ¼ ð1; þ1Þ, we use Hermite polynomials with 2 wðsÞ / es =2 . Upon substituting (4.4) into (4.3) and enforcing the residue to be orthogonal to the linear polynomial space of degree up to n, we obtain the following Galerkin system n X ^j du ^i ; ¼ eij u dt i¼0
^ j ð0Þ ¼ d0;j ; u
ð4:6Þ
R where eij ¼ IX sUi ðsÞUj ðsÞwðsÞds. These are a set of ODEs that can be easily solved. And the resulting numerical solution un converges to u in L2w norm. (See Appendix for proof.) We then assign various distributions to Z a posterior, use these distributions to approximate the statistics of the solution via un directly (without resorting to the numerical simulations), and examine the errors induced by this procedure.
4.1.1. Bounded encapsulation We first assume the exact (and unknown) distribution of Z is an exponential distribution, i.e., qZ ðsÞ expðsÞ. For the encapsulation problem (4.3) we use X 2 ð0; bÞ with various length of b. This is the similar case considered in [9], where the convergence in mean, a weak error measure, was discussed. Here we study the convergence in L2q norm, a strong norm, and examine its behavior based on Theorem 3.2. The results are depicted in Fig. 4.1. It is obvious that the errors converge (exponentially fast) as the order of expansions is increased, before they saturate at certain level. The error saturation level depends on the value of b — a larger value of b induces a smaller error because it causes a smaller ‘‘truncation’’ in the distribution of Z and hence a smaller value d in the second error term in (3.23) of Theorem 3.2. Weak error measure, e.g., error in mean value, behaves similarly. It has been documented in [9] and is not presented here. We further examine the error dependence on b, which is a direct indication of the ‘‘truncation’’ error d. In Fig. 4.2, the strong and weak errors are examined against the parameter b, at a fixed and sufficiently high polynomial order of n ¼ 18. Numerical errors in the Galerkin solver are thus negligible. It can be seen that the errors decrease when b increases, as they should. It is also evident that the strong error converges slower than the weak error with respect to b. This is consistent with the estimate from (3.19) and (3.23), as the strong error scales as d1=p (p ¼ 2 in this case) but weak error scales as d.
218
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
Fig. 4.1. Convergence of the errors in strong L2q norm for linear ODE (4.1) at t ¼ 1, where the encapsulation set is X 2 ð0; bÞ and the exact distribution is exponential qZ ðsÞ expðsÞ.
Fig. 4.2. Convergence of the errors of the strong norm and the weak norm for linear ODE (4.1) at t ¼ 1, where the encapsulation set is X 2 ð0; bÞ and the exact distribution is exponential qZ expðsÞ. The results are for fixed polynomial order of n ¼ 18 and various values of b.
4.1.2. Semi-bounded encapsulation We now employ X 2 ð0; þ1Þ as the encapsulation variable. The corresponds to the practical case where one does not have a reliable estimate of the upper bound of the epistemic variable. We then use Laguerre polynomial with weight wðsÞ / es to approximate the encapsulation problem (4.3). We first assume the true distribution of Z is uniform in (0,1) with qZ ðsÞ ¼ 1. In this case there is no ‘‘truncation’’ error because IX ¼ ð0; þ1Þ completely encapsulates the real variable Z. Error convergence, in both the strong form of L2q and the weak form in mean, are plotted in Fig. 4.3, for different orders of polynomial approximations at a fixed time t ¼ 1. It can be seen that in this case both the strong error and the weak error converge at the same rate – exponential rate – as
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
219
Fig. 4.3. Convergence of errors in the strong form of L2q norm and the weak form of mean, along with the L2w error (n ) of the Galerkin solutions of (4.3), at different polynomial approximation orders n. Here X 2 ð0; 1Þ and the true a posterior distribution of Z is uniform in (0, 1).
the L2w polynomial approximation error of (4.3). This is consistent with the error bounds from Theorems 3.1 and 3.2, as in this case the complete encapsulation of Z by X ensures d ¼ 0. We now assume the true distribution of Z is uniform, but resides in different intervals of (a,b). In Fig. 4.4, the errors are presented from intervals of (0, 1), (5, 6), and (10, 11), in both the L2q form (left) and mean (right). It is obvious that all errors converge (exponentially fast). It also should be remarked that although all intervals have the same size of one, as the interval becomes further away from the origin the errors become larger. This is because in the Laguerre approximation the weight is wðsÞ / es — it achieves the best accuracy close to zero and becomes progressively less accurate away from zero. Consequently, we should emphasize that even though the ability to use semi-bounded and unbounded encapsulation variables in the current framework does alleviate the modeling effort, sometimes significantly, to estimate the ranges of the epistemic variables, it does not imply this can completely eliminate the need to estimate the epistemic variables. When using unbounded domains, the accuracy of the numerical approximation will inevitably be less satisfactory in the region where the weight function approaches zero. It is thus important to have a good, or even just a vague, idea of the true range of the epistemic variables so that one can put the estimated range close to the ‘‘peak’’ of the weight function to ensure more accurate solutions.
Fig. 4.4. Convergence of errors for uniformly distribution Z in different intervals of (a, b), using encapsulation variable X 2 ð0; þ1Þ. Left: L2q error; Right: error in mean.
220
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
4.1.3. Unbounded encapsulation Finally we examine a case where the encapsulation variable is assumed to be X 2 ð1; 1Þ. Correspondingly we employ Hermite polynomial to solve (4.3) numerically. This is a simplified case of the practical situation where one does not have a good understanding of either the upper bound or the lower bound of the epistemic variables. Using an unbounded encapsulation variable would serve as one way to circumvent the difficulty. Both the errors in the strong form (L2q ) and the weak form (mean) are shown in Fig. 4.5, along with the L2w convergence error of the numerical solution of (4.3). Once again we emphasize that in practice it is desirable that one has at least a vague idea of where the true range of the epistemic variable is, and can then ‘‘center’’ the weight function wðsÞ of the Hermite polynomial approximation around that range to achieve more accurate numerical results. 4.2. Homogenous Diffusion Equation We consider the following stochastic diffusion problem in one spatial dimension and multiple random dimensions (d > 1).
@ @x
jðx; ZÞ
@u ðx; ZÞ ¼ f ðx; ZÞ; @x
ðx; ZÞ 2 ð0; 1Þ Rd :
ð4:7Þ
Without loss of generality we impose simple boundary conditions
uð0; ZÞ ¼ 0;
uð1; ZÞ ¼ 0;
Z 2 Rd
and let f ¼ 2 be a constant. For the diffusivity field, we employ the following construction d X 1
jðx; ZÞ ¼ 1 þ r
k¼1
k
2
p2
ð1 þ cosð2pkxÞÞZ k ;
ð4:8Þ
where r > 0 is a parameter to control the variation level of the diffusivity field. Here we choose this fixed form representation to focus on the treatment of the random variables Z k ’s, which are assumed to be epistemic. Note that usually the expansion is obtained by certain decomposition procedure, e.g., Karhunen–Loeve expansion. For non-Gaussian processes, the decomposition usually can not fully determine the probability distributions of Z k ’s. Since the detailed examination of using different encapsulation strategies has been carried out in the previous example, here we mostly focus on dealing with the unknown dependence structure of the epistemic random variables Z k ’s. We set d ¼ 6 and r ¼ 1 in the expansion (4.8). (These choices are rather arbitrary.) We assume that the true (and unknown) distributions of the epistemic variables are: Z 1 betað0; 1; 3; 2Þ, and Z 3 betað1; 0; 1; 1Þ, Z 5 betað0:5; 0:5; 0; 0Þ, and Z 2 ¼ Z 1 Z 5 , Z 4 ¼ ðZ 21 þ 1ÞZ 3 , and Z 6 ¼ Z 5 . Clearly all the variables are bounded, Z 1 2 ð0; 1Þ;
Fig. 4.5. Convergence of the errors of in the strong form of L2q norm and the weak form of mean, along with the L2w error (n ) of the Galerkin solutions of (4.3), at different polynomial approximation orders n. Here the encapsulation set is X 2 ð1; 1Þ and the exact distribution is uniform in (0, 1).
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
221
Fig. 4.6. Convergence of the errors of in both strong form (L2 norm) and in a weak form (in mean) for the diffusion problem. Errors are computed against the high-order (n ¼ 8) Galerkin solution at a fixed spatial location x ¼ 0:5257.
Z 2 2 ð0:5; 0:5Þ, Z 3 2 ð1; 0Þ; Z 4 2 ð2; 0Þ; Z 5 2 ð0:5; 0:5Þ, and Z 6 2 ð0:5; 0:5Þ. Also by construction the variables are dependent. What appears to be a six-dimensional problem is in fact a three-dimensional one, for Z 2 ; Z 4 and Z 6 are functions of the other three variables. In the following epistemic analysis we assume that none of the above information regarding the distributions is unavailable. The only information we have is a conservative estimate of the bounds of the variables. And based on the bound estimation, we employ a straightforward linear transformation and define X 2 IX ¼ ð0; þ1Þ6 to completely encapsulate the epistemic variables Z. Our encapsulation problem is therefore
@ @u ½jðx; XÞ ðx; XÞ ¼ f ðx; XÞ; @x @x
ðx; XÞ 2 ð0; 1Þ ð0; þ1Þ6 ;
ð4:9Þ
where same boundary conditions and source term f hold, and the diffusivity field follows the similar construction of (4.8) with the encapsulation variables X in place of Z. Clearly the encapsulation problem is six-dimensional in the parameter space of X. We employ the six-dimensional Laguerre polynomials, in conjunction with a Galerkin procedure, with weight function Q wðsÞ / 6k¼1 esk to solve the problem. Due to lack of analytical solution, here we employ a sufficiently high-order polynomials solution at n ¼ 8 order as the ‘‘numerical exact’’ solution, against which we compare the solutions obtained by lower-order expansions. Both the strong errors (in L2q ) and weak errors (in mean value) are computed and presented in Fig. 4.6. The errors are computed at a fixed spatial location (x ¼ 0:5257), where the solution is non-trivial. (One can of course take a normed error in the spatial dimension as well.) The usual (exponentially) fast convergence with respect to the approximation order can be clearly observed. No error saturation is present because the encapsulation variables can fully encapsulate the epistemic variables. Various other tests have also be carried out to examine different encapsulation strategies, especially when the encapsulation variables ‘‘truncate’’ the true ranges of the epistemic variables. The results are quite similar to those presented in the previous example, and have also been detailed in the work of [9]. Hence they are not presented here.
5. Summary In this paper we present a computational strategy for epistemic uncertainty analysis. The new method is a notable extension of the work of [9], in the sense that it allows one to use unbounded intervals to encapsulate the true (and unknown) ranges of the epistemic variables. Consequently the new method can employ numerical solutions that are accurate in Lp norms. These feature makes the method more flexible in practice. Theoretical analysis of the method is conducted, where errors in both strong form and weak form are analyzed. Numerical tests are conducted to verify the theory. Though more extensive work is required, and is ongoing, to further examine the methodology, the new method appears to be a fairly effective tool for epistemic uncertainty analysis.
222
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
Appendix: Convergence of Galerkin solution for ODE We now present the convergence proof for the orthogonal polynomial Galerkin solution for the ordinary differential Eq. (4.3). In fact, we consider a slightly more general form of the problem
du ¼ aðXÞu; dt
0 < t 6 T;
uðt ¼ 0Þ ¼ 1;
ð5:1Þ
where the coefficient a is bounded away from infinity
1 < amin 6 a 6 amax < þ1:
ð5:2Þ
Let IX be the range of X and fUj g be a set of orthogonal polynomials defined on IX
Z
Uj ðsÞUk ðsÞwðsÞds ¼ djk ;
ð5:3Þ
IX
R where w > 0 is a weight function inducing an inner product ðf ; gÞw ¼ f ðsÞgðsÞwðsÞds and norm kf kw ¼ ðf ; f Þ1=2 w . Here and m after, the space Hw ðIX Þ denotes the usual weighted Sobolev spaces equipped with the norm kf kHmw . In practice, it is often convenient to map IX to a standard interval such as ð1; þ1Þ, ð0; þ1Þ, or ð1; 1Þ, where standard forms of orthogonal polynomial can be readily applied. Hereafter we will not explicitly specify the interval. In Gelerkin method we seek an expansion
v N ðt; XÞ ¼
N X
v^ j ðtÞUj ðXÞ
ð5:4Þ
j¼0
such that (5.1) is satisfied in a weak form, i.e., the residue is orthogonal to the polynomial space where the solution is constructed. This leads to, for all k ¼ 0; . . . ; N,
dv N ¼ 0: þ av N ; Uk dt w Straightforward derivation then gives us, for all k ¼ 0; . . . ; N, N X dv^ k ¼ ajk v^ j ; dt j¼0
ð5:5Þ
where
ajk ¼
Z
aðsÞUj ðsÞUk ðsÞwðsÞds
ð5:6Þ
^ k ðt ¼ 0Þ ¼ dk0 . and the initial conditions are v In order to prove the convergence of the Galerkin solution v N to the solution u of (5.1), it is convenient to consider an intermediate projection. More precisely, for uðt; XÞ 2 L2w ðIX Þ for any t, its orthogonal projection P N u exists and is defined as
PN uðt; XÞ ,
N X
^ j ðtÞUj ðXÞ; u
^ j ¼ ðu; Uj Þw : u
ð5:7Þ
j¼0
Obviously, P N u is optimal in the k kw norm and limN!1 P N u ¼ u. Assume that for some m > 1=2,
kðPN u uÞðtÞkw 6 CN m kuðtÞkHmw ðIX Þ
8 t > 0;
ð5:8Þ
where the constant C does not depend on N. Note that the rate of convergence depends on the regularity of u and the type of orthogonal polynomials fUj g. Theorem 5.1 (Convergence). Suppose that the assumption (5.8) holds. The Galerkin solution v N from (5.4)–(5.6) converges to the solution u of the Eq. (5.1). Moreover, if the solution uðt; XÞ belongs to the weighted Sobolev space Hm w ðI X Þ for any t, there exists a constant C, independent of N, such that
maxkðv N uÞðtÞkw 6 CNmþ1=2 maxkuðtÞkHmw ðIX Þ :
06t6T
06t6T
ð5:9Þ
^ j g given by (5.7) satisfy the following equations, Proof. It is easy to verify that the coefficients fu 1 X ^k du ^j ; ¼ ajk u dt j¼0
k ¼ 0; 1; . . . ; 1:
ð5:10Þ
223
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
Let us define
^k ; ^ek ¼ v^ k u
k ¼ 0; . . . ; N:
ð5:11Þ
We then have N 1 X X d^ek ^j : ajk u ¼ ajk ^ej þ dt j¼0 j>N
Upon multiplying both sides by ^ ek and taking a summation over k, we obtain N N X N N X 1 X X X d^e ^ek k ¼ ^ek ajk ^ej þ ^ek ajk u ^j : dt k¼0 k¼0 j¼0 k¼0 j>N
Let eN ðt; XÞ ¼ v N P N u. Naturally we have keN kw ¼
e ¼ ð^e0 ; . . . ; ^eN ÞT ;
PN
^2 k¼0 ek .
We now adopt vector–matrix notation by defining
A ¼ ðajk Þ06j;k6N :
The above equation can then be written as
1 d kek2 ¼ eT Ae þ eT r; 2 dt
ð5:12Þ
where the norm is the vector 2-norm and kek ¼ keN kw and r ¼ ð^r0 ; . . . ; ^r N ÞT with the entries
^rk ¼
1 X ^j ; ajk u
k ¼ 0; . . . ; N:
ð5:13Þ
j>N
Lemma 5.2. Following the definitions of e and A and the condition (5.2)
eT Ae P amin kek2 :
ð5:14Þ
Proof. Consider an arbitrary vector e and construct a function
qðsÞ ¼ ^e0 þ ^e1 U1 ðsÞ þ þ ^eN UN ðsÞ: Obviously qðsÞ is a polynomial and satisfies kqkw ¼ kek.
eT Ae ¼
Z N X N N X N Z X X ^ej ajk ^ek ¼ aðsÞ^ej ^ek Uj ðsÞUk ðsÞwðsÞds ¼ aðsÞq2 ðsÞwðsÞds P amin kqk2w : j¼0 k¼0
ð5:15Þ
j¼0 k¼0
h Substituting this result into (5.12), we obtain
d kek 6 amin kek þ krk: dt
ð5:16Þ
We now analyze the term krk. Using the definition (5.13), we derive
^rk ¼
Z Z 1 1 Z 1 X X X ^j ¼ ajk u aðsÞUj ðsÞUk ðsÞwðsÞdsu^ j ¼ aðsÞ u^ j Uj ðsÞUk ðsÞwðsÞds ¼ aðsÞðu PN uÞUk ðsÞwðsÞds: j>N
j>N
j>N
Then, it follows that
krk ¼
N X ^r2k k¼0
!1=2 ¼
N Z X k¼0
2 !1=2 6 ðu PN uÞaðsÞUk ðsÞwðsÞds
pffiffiffiffiffiffiffiffiffiffiffiffiffi 6 amax N þ 1ku PN ukw ; , bðN; tÞ;
N X ku PN uk2w kaUk k2w
!1=2
k¼0
ð5:17Þ
where the Cauchy–Schwarz inequality and the orthonormal property of the basis are used. Note that for m > 1=2; bðN; tÞ ! 0 as N ! 1 for all t due to (5.8). And from (5.16 and 5.17), we obtain
d kek 6 amin kek þ bðN; tÞ: dt
ð5:18Þ
224
X. Chen et al. / Journal of Computational Physics 240 (2013) 211–224
The Gronwall inequality then leads to
keðtÞk 6 eamin t keð0Þk þ
Z
t
eamin ðtsÞ bðN; sÞds:
ð5:19Þ
0
An application of the triangle inequality implies
kv N ukw 6 kv N PN ukw þ kPN u ukw and the convergence theorem is then established for N ! 1 by using (5.8), (5.17) and keð0Þk ¼ 0.
h
It is trivial to see that the same convergence result holds true for the case of keð0Þk – 0 but converges at the same rate of the projection error ku PN ukw . This corresponds to the case of (5.1) with a random initial condition u0 ðXÞ that is approximated by its orthogonal projection PN u0 . Correspondingly the initial conditions of the system (5.5) becomes v^ k ðt ¼ 0Þ ¼ ðu0 ; Uk Þw . This is the typical procedure in practice. References [1] I. Babuska, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal. 45 (3) (2007) 1005–1034. [2] I. Babuska, R. Tempone, G.E. Zouraris, Galerkin finite element approximations of stochastic elliptic differential equations, SIAM J. Numer. Anal. 42 (2004) 800–825. [3] D. Dubois, H. Prade, Possibility Theory: An Approach to Computerized Processing of Uncertainty, Plenum, New York, 1998. [4] M. Eldred, Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design, in: 50th AIAA/ ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, vol. AIAA-2009-2249, 2009. [5] M. Eldred, L. Swilder, G. Tang, Mixed aleatory-epistemic uncertainty quantification with stochastic expansions and optimization-based interval estimation, Reliab. Eng. Syst. Saf. 96 (9) (2011) 1092–1113. [6] R.G. Ghanem, P. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer-Verlag, 1991. [7] J.C. Helton, J.D. Johnson, W.L. Oberkampf, C.J. Sallaberry, Representation of analysis results involving aleatory and epistemic uncertainty, Technical Report 4379, Sandia National Laboratories, 2008. [8] J.C. Helton, J.D. Johnson, W.L. Oberkampf, C.B. Storlie, A sampling-based computational strategy for the representation of epistemic uncertainty in model predictions with evidence theory, Technical Report 5557, Sandia National Laboratories, 2006. [9] J. Jakeman, M. Eldred, D. Xiu, Numerical approach for quantification of epistemic uncertainty, J. Comput. Phys. 229 (2010) 4648–4663. [10] L. Swiler, T. Paez, R. Mayes, M. Eldred. Epistemic uncertainty in the calculation of margins, in: AIAA Structures, Structural Dynamics, and Materials Conference, Palm Springs CA, May 2009. [11] L.P. Swiler, T.L. Paez, R.L. Mayes, Epistemic uncertainty quantification tutorial, in: IMAC XXVII Conference and Exposition on Structural Dynamics, Orlando, FL, Feburary 2009. [12] D. Xiu, Efficient collocational approach for parametric uncertainty analysis, Commun. Comput. Phys. 2 (2) (2007) 293–309. [13] D. Xiu, Fast numerical methods for stochastic computations: a review, Commun. Comput. Phys. 5 (2009) 242–272. [14] D. Xiu, Numerical Methods for Stochastic Computations, Princeton Univeristy Press, Princeton, New Jersey, 2010. [15] D. Xiu, J.S. Hesthaven, High-order collocation methods for differential equations with random inputs, SIAM J. Sci. Comput. 27 (3) (2005) 1118–1139. [16] D. Xiu, G.E. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput. 24 (2) (2002) 619–644.