51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference
18th 12 - 15 April 2010, Orlando, Florida
AIAA 2010-2924
Simplex Elements Stochastic Collocation in Higher-Dimensional Probability Spaces Jeroen A.S. Witteveen∗, Gianluca Iaccarino† Center for Turbulence Research, Stanford University, Building 500, Stanford, CA 94305–3035, USA
A Simplex Elements Stochastic Collocation (SESC) method is introduced for robust and efficient propagation of uncertainty through computational models. The presented non– intrusive Uncertainty Quantification (UQ) method is based on adaptive grid refinement of a simplex elements discretization in probability space. The approach is equally robust as Monte Carlo (MC) simulation in terms of the Extremum Diminishing (ED) robustness concept. The efficiency of SESC is based on high degree polynomial interpolation, randomized refinement sampling, and Essentially Extremum Diminishing (EED) extrapolation. This results in a superlinear convergence rate and a linear increase of the initial number of samples with dimensionality. The flexibility of simplex elements is further employed to discretize non–hypercube probability spaces with correlated random parameters.
I.
Introduction
tate-of-the-art uncertainty quantification methods are based on a multi-element discretization of probS ability space. The piecewise polynomial approximation of these multi-element methods is robust for a large class of problems. This is an essential property for using uncertainty quantification to improve 3–5, 11–13, 17
the confidence in numerical predictions. Multi–element discretizations are significantly more reliable than Stochastic Collocation1, 6, 15, 25 methods based on global polynomial Lagrangian interpolation through samples in selected Gauss quadrature points in probability space. This non–intrusive version of the Stochastic Galerkin7, 18, 24 method has been developed as a more efficient alternative for the optimally robust Monte Carlo8 simulation method owing to its spectral convergence for smooth responses. However, it is widely acknowledged14, 26 that for discontinuities or numerical noise in the response surface, Stochastic Collocation can result in an oscillatory approximation due to its global polynomial interpolation. On the other hand, multi–element approaches are usually also based on employing single–element Stochastic Collocation independently in multiple hypercube elements discretizing the probability space. For higher degree interpolations these methods can still result in local oscillations in elements that contain a discontinuity, and overshoots can be present even for low degree approximations. Often not all samples in an element can be reused after refinement and tensor product extensions to higher dimensions may be employed, which compromises the efficiency of multi–element discretizations. In this paper, a Simplex Elements Stochastic Collocation (SESC) method is introduced that combines the robustness of Monte Carlo (MC) simulation with the superlinear convergence of Stochastic Collocation (SC). SESC is a significant extension of the Adaptive Stochastic Finite Elements (ASFE) method based on Newton–Cotes quadrature in simplex elements19, 20 from first and second degree quadrature to higher degree polynomial interpolation, randomized refinement sampling, and Essentially Extremum Diminishing (EED) extrapolation. These developments make SESC a fundamentally multi–dimensional approach that is more efficient than ASFE for multiple uncertainties. Simplex elements are employed here based on the viewpoint that optimal uncertainty quantification (UQ) methods for single–element applications are not necessarily the ∗ Postdoctoral † Assistant
Fellow, Member AIAA, Phone: +1 650 723 9601, Fax: +1 650 723 9617,
[email protected]. Professor, Member AIAA,
[email protected].
1 of 17 American Institute of Aeronautics Copyright © 2010 by J.A.S. Witteveen, G. Iaccarino. Published by the American Institute of Aeronauticsand and Astronautics Astronautics, Inc., with permission.
optimal choice for multi–element discretizations. Simplices have in a multi–element formulation the advantage that they lead to a UQ method that is equally robust as MC in terms of the Extremum Diminishing10, 21 (ED) robustness concept in probability space. The combination with Newton–Cotes quadrature also limits the number of samples owing to the location of a considerable number of the Newton–Cotes points on the boundaries of simplex elements, such that samples are shared by adjacent elements. The introduction of randomized refinement sampling into the SESC method instead of higher degree Newton–Cotes quadrature results in a superlinear convergence behavior. To that end, a simplex is refined by randomly selecting a new sampling point in the element instead of using a deterministic refinement rule. This employs the flexibility of a simplex elements discretization to interpolate scattered sampling points using a Delaunay triangulation. The number of samples in the initial discretization is also reduced to a linear increase with dimension by including extrapolation in the SESC formulation. This leads to the introduction of a quantitative robustness measure in terms of the Essentially Extremum Diminishing (EED) concept, which is derived from the Local Extremum Diminishing (LED) property extended to probability space in this paper. SESC further uses a solution–based refinement criterion to adaptively refine one element at a time, where all samples are reused after the hierarchical refinements. The refinement is automatically terminated when an error estimate reaches a stopping criterion. The flexibility of simplex elements is also employed to discretize non–hypercube probability spaces. The general formulation of SESC is outlined in section II. The randomized refinement sampling formulation based on Delaunay triangulation is developed in section III. Error estimates are derived and tested in a numerical application to an analytical test function. In section IV the method is extended to EED extrapolation and applied to the uniform distribution, the non–uniform beta distribution, and the normal distribution with an unbounded support. The flexibility of the simplex elements discretization with extrapolation is also employed to treat non–hypercube probability spaces with correlated random parameters. Conclusions are summarized in section V.
II.
Simplex Elements Stochastic Collocation
A general formulation of Simplex Elements Stochastic Collocation (SESC) is considered as a reference for the extensions in subsequent sections. A.
General formulation
Consider the following computational problem for an output of interest u(x, t, ξ(ω)) L(x, t, ξ(ω); u(x, t, ξ(ω))) = S(x, t, ξ(ω)),
(1)
with appropriate initial and boundary conditions. The operator L and source term S are defined on domain D × T × Ξ, where x ∈ D and t ∈ T are the spatial and temporal dimensions with D ⊂ Rd , d ∈ {1, 2, 3}, and T ⊂ R. Randomness is introduced in (1) and its initial and boundary conditions in terms of nξ second– order random parameters ξ(ω) = {ξ1 (ω1 ), . . . , ξnξ (ωnξ )} ∈ Ξ with parameter space Ξ ⊂ Rnξ . The symbol ω = {ω1 , . . . , ωnξ } ∈ Ω ⊂ Rnξ denotes realizations in the probability space (Ω, F, P ) with F ⊂ 2Ω the σ– algebra of events and P a probability measure. The random variables ω are by definition standard uniformly distributed as U(0, 1). Random parameters ξ(ω) can have any arbitrary probability density fξ (ξ(ω)). For parameter space Ξ then holds fξ (ξ(ω)) > 0 for all ξ(ω) ∈ Ξ and fξ (ξ(ω)) = 0 for ξ(ω) ∈ / Ξ. The argument ω will be dropped from here on to simplify the notation. The aim of uncertainty propagation is then to find the probability distribution of u(x, t, ξ) and its statistical moments µui (x, t) given by Z u(x, t, ξ)i fξ (ξ)dξ. (2) µui (x, t) = Ξ
A multi–element UQ method computes this weighted integral as a summation of integrals over ne disjoint elements Ξ = Ξ1 ∪ · · · ∪ Ξne ne Z X µui (x, t) = u(x, t, ξ)i fξ (ξ)dξ. (3) j=1
Ξj
In SESC the integrals in the elements are computed by approximating response surface u(ξ) by an interpolation w(ξ) of ns samples v = {v1 , . . . , vns }. Here the arguments x and t are omitted for clarity of the 2 of 17 American Institute of Aeronautics and Astronautics
notation. Non–intrusive SESC uncertainty quantification method q then consists of a sampling method g and an interpolation method h, for which holds w(ξ) = q(u(ξ)) = h(g(u(ξ))). The sampling method g selects the sampling points ξ k for k = 1, . . . , ns and returns the sampled values v = g(u(ξ)), with vk = gk (u(ξ)) = u(ξ k ). Sample vk is computed by solving (1) for realization ξ k of the random parameter vector ξ L(x, t, ξ k ; vk (x, t)) = S(x, t, ξ k ),
(4)
for k = 1, . . . , ns . The interpolation of the samples w(ξ) = h(v) consists of a piecewise polynomial function for ξ ∈ Ξj ,
w(ξ) = wj (ξ),
(5)
with wj (ξ) the polynomial interpolation of degree p of the samples vj = {vkj,0 , . . . , vkj,N } at the sampling points {ξkj,0 , . . . , ξ kj,N } in element Ξj , where kj,l ∈ {1, . . . , ns } for j = 1, . . . , ne and l = 0, . . . , N with N +1=
(nξ + p)! , nξ !p!
(6)
the number of samples in the elements. The polynomial interpolation wj (ξ) in element Ξj can then be expressed in terms of a truncated Polynomial Chaos expansion7, 18 wj (ξ) =
a0 Γj,0 +
nξ X
ai1 Γj,1 (ξi1 ) +
nξ i1 X X
i1 =1 i2 =1
ai1 ,i2 Γj,2 (ξi1 , ξi2 ) +
i1 =1 i2 =1
i1 =1
...+
nξ i1 X X
ip−1
···
X
ai1 ,...,ip Γj,p (ξi1 , . . . , ξip ),
(7)
ip =1
with multi–dimensional basis polynomials Γj,p˜ of exact degree p˜. Expansion (7) can be written in the following more convenient shorter notation wj (ξ) =
N X
cj,l Ψj,l (ξ),
(8)
l=0
with a one–to–one correspondence between the basis polynomials Γj,p˜ and Ψj,l , and the coefficients ai1 ,...,ip˜ and cj,l . If the response approximation wj (x, t, ξ) depends in addition to the random parameters ξ also on spatial x and time t coordinates, then there occurs a separation of variables in terms of the random parameters, Ψj,l (ξ), and the spatial and temporal dimensions, cj,l (x, t). The polynomial coefficients cj,l can be determined from the interpolation condition wj (ξ kj,l ) = vkj,l , for l = 0, . . . , N , which leads to Ψj,0 (ξ kj,0 ) Ψj,0 (ξ kj,1 ) .. . Ψj,0 (ξ kj,N )
(9)
the matrix equation9, 16 Ψj,1 (ξ kj,0 ) · · · Ψj,N (ξ kj,0 ) Ψj,1 (ξ kj,1 ) · · · Ψj,N (ξ kj,1 ) .. .. .. . . . Ψj,1 (ξ kj,N ) · · · Ψj,N (ξ kj,N )
cj,0 cj,1 .. . cj,N
=
vkj,0 vkj,1 .. . vkj,N
.
(10)
The piecewise polynomial approximation w(ξ) of response surface u(ξ) is eventually found by substituting (8) into (5) with cj,l from (10) for j = 1, . . . , ne and l = 0, . . . , N . The probability distribution function and the statistical moments µui of u(ξ) given by (3) are then approximated by the probability distribution and the moments µwi of w(ξ) ne Z X µui (x, t) ≈ µwi (x, t) = wj (x, t, ξ)i fξ (ξ)dξ, (11) j=1
Ξj
in which the multi–dimensional integrals are evaluated using a weighted Monte Carlo integration of response surface approximation w(ξ) with nsMC ≫ ns integration points. This is a fast operation, since it only involves sampling of piecewise polynomial function w(ξ) given by (5) and no additional evaluations of the exact response u(ξ). We have not encountered difficulties in solving (10) using Matlab, although the matrix can have a high condition number. 3 of 17 American Institute of Aeronautics and Astronautics
III.
Randomized refinement sampling
Random sampling at the refinement of an element is introduced in the general SESC formulation to obtain a superlinear convergence behavior. This is motivated by the knowledge that random sampling is an effective sampling strategy that naturally avoids a tensor grid of samples in multiple dimensions. The flexibility of the simplex elements discretization is employed to interpolate the randomized samples using a Delaunay triangulation. Higher degree polynomial interpolation is obtained using a stencil of samples of surrounding elements. Error estimates for this formulation are assessed in application to an analytical test function. A.
Delaunay triangulation
The Delaunay triangulation2 of the random samples maximizes the minimum angles in the simplex elements discretization of parameter space. It is used as the basis for interpolating the scattered sampling points with a higher degree polynomial. 1.
Randomized element refinement
The samples ξk in the initial discretization are located in the vertices of the hypercube parameter space Ξ and one in the interior only. The elements are subsequently refined by sampling a random location ξk ∈ Ξj in the element. In order to ensure a sufficient spread of the samples, the random sample is confined to a simplex sub–element ξ k ∈ Ξsubj ⊂ Ξj . The nξ + 1 vertices ξ subj,l of the nξ –simplex Ξsubj , with l = 0, . . . , nξ , are defined as the centers of the (nξ − 1)–faces of simplex Ξj ξsubj,l
nξ 1 X = ξ kj,l∗ , nξ ∗
(12)
l ∗=0 l 6=l
with kj,l∗ ∈ {1, . . . , ns } and j = 1, . . . , ne , see Figure 1a for a two–dimensional example. In one dimension nξ = 1, this reduces to sampling in the center of element Ξj . The Delaunay triangulation of the samples is reconstructed after each element refinement while reusing all previous samples. This leads to a discretization that combines randomness with a good spread of the samples shown in Figure 1b for ns = 17. The measure used for selecting an element Ξj for refinement in this approach is discussed in section IV.C on non–uniform distributions.
0.8
0.8
0.6
0.6
ξ2
1
ξ2
1
0.4
0.4
0.2
0.2
0 0
0.2
0.4
ξ
0.6
0.8
0 0
1
0.2
0.4
ξ
0.6
0.8
1
1
1
(a) Initial refinement with ns = 6
(b) Refined mesh with ns = 17
Figure 1. Discretization of two–dimensional parameter space Ξ with nξ = 2 using randomized refinement sampling with simplex sub–element Ξsubj given by the dashed lines.
4 of 17 American Institute of Aeronautics and Astronautics
2.
Higher degree polynomial interpolation
Higher degree polynomial interpolation is essential for efficient UQ of multiple uncertainties. Therefore, the polynomial interpolation wj (ξ) (8) N X wj (ξ) = cj,l Ψj,l (ξ), l=0
in element Ξj , is constructed using a stencil Sj consisting of the samples in the vertices of surrounding elements. The stencil of N + 1 sampling points Sj = {ξkj,0 , . . . , ξ kj,N } and the corresponding samples {vkj,0 , . . . , vkj,N }, with kj,l ∈ {1, . . . , ns }, is then used in (10) to find the polynomial coefficients cj,l for j = 1, . . . , ne and l = 0, . . . , N with N given by (6) N +1=
(nξ + p)! . nξ !p!
The stencil Sj is constructed as follows based on the nearest neighbor principle. The first nξ + 1 points {ξkj,0 , . . . , ξkj,n } are the vertices of element Ξj , see Figure 2a. This stencil of N +1 = nξ +1 sampling points ξ corresponds to a piecewise linear interpolation of degree p = 1. For higher degree polynomial interpolation, surrounding sampling points ξk are added to Sj in the order of their Euclidean distance to the center ξ centerj of simplex Ξj (13) ||ξ k − ξ centerj ||, with k ∈ {1, . . . , ns }\{kj,0 , . . . , kj,nξ }, where the center of element Ξj with vertices ξ kj,l for l = 0, . . . , nξ is defined as nξ 1 X ξcenterj = ξ kj,l . (14) nξ + 1 l=0
Stencils for p = 2 and p = 4 with N + 1 = 6 and N + 1 = 15 are shown in Figures 2b and 2c. The polynomial degree p of interpolation wj (ξ) is in this formulation limited by the total number of samples ns following (6) (nξ + p)! = N + 1 ≤ ns . nξ !p!
(15)
Degree p can then be increased with increasing ns such that it satisfies (15) and (nξ + (p + 1))! > ns . nξ !(p + 1)!
(16)
It can also be limited to a maximum degree p ≤ pmax .
(a) p = 1 with N + 1 = 3
(b) p = 2 with N + 1 = 6
(c) p = 4 with N + 1 = 15
Figure 2. Higher degree SESC stencils Sj in two–dimensional parameter space Ξ with nξ = 2.
The robustness of the SESC method based on randomized refinement sampling is guaranteed by the Local Extremum Diminishing23 (LED) and Extremum Diminishing21 (ED) properties extended to probability space. The first degree formulation with p = 1 is by definition LED in all elements Ξj for j = 1, . . . , ne . Since the elements Ξj cover the whole hypercube probability space Ξ, the method is also ED. For higher
5 of 17 American Institute of Aeronautics and Astronautics
degree interpolations it is verified whether wj (ξ) is LEC with respect to the samples vj in the vertices of element Ξj . If the Monte Carlo integration for computing the statistical moments violates min wj (ξ MCk ) ≥ min vj ∧ max wj (ξ MCk ) ≤ max vj , l
l
j,l
j,l
(17)
for l = 1, . . . , nsMCj , then the local polynomial degree pj of wj (ξ) in Ξj and the corresponding stencil Sj are reduced until (17) is satisfied with (nξ + pj )! . (18) Nj + 1 = nξ !pj ! This check is repeated for all elements involved in stencil Sj and for all elements Ξj . LEC limiter (17) is then automatically used both for robustly approximating discontinuities and for suppressing potential Runge phenomena in the higher degree polynomial interpolation of the scattered sampling points. B.
Error estimates
Error estimates for SESC with randomized refinement sampling are based on expression22 ε(ξ kj,ref ) = wj (ξ kj,ref ) − vkj,ref ,
(19)
for the error in approximation wj (ξ) in the new sampling point ξ kj,ref at the refinement of element Ξj . Sampling point ξkj,ref is a random location in sub–simplex Ξsubj . This leads in each sampling point ξ k to a known error εk = w(ξ k ) − vk between sample vk and the approximation w(ξ) prior to adding the sampling point at ξ k . The approximation of the error ε˜j in element Ξj is then given by the error εk at the most recently added sampling point ξ k out of the vertices {ξkj,0 , . . . , ξkj,n } of the simplex ξ
ε˜j = εl ,
(20)
l = max{kj,0 , . . . , kj,nξ }.
(21)
with The relation between approximation ε˜j and error estimate εˆj after refinement is ε˜j
εˆj = 2
pj +1 nξ
,
(22)
where pj is the current polynomial degree of approximation wj (ξ) in element Ξj . This formulation accounts for a potentially different pj in different elements due to the LEC constraint (17), and for the increasing degree pj with increasing number of samples ns . This leads to the following error estimates for SESC with randomized refinement sampling for the mean εˆµ X ne ε˜j 1 n (23) εˆµ = s , nsMC j=1 MCj pnj +1 ξ 2 and for standard deviation εˆσ εˆσ
=
nsMC j ne X X 1 σw − nsMC j=1 ˜ k=1 nsMC
ne Xj X 1 nsMC j=1 ˜ k=1
and L2 error norm εˆL2 εˆL2
w(ξ MCk
˜ j,k
ε˜j
)− 2
pj +1 nξ
!2
−
1 !2 2 ε˜j w(ξ MCk ) − pj +1 , ˜ j,k 2 nξ
v u ne u 1 X ns =t nsMC j=1 MCj
ε˜j 2
pj +1 nξ
!2
.
6 of 17 American Institute of Aeronautics and Astronautics
(24)
(25)
C.
Superlinear convergence
Results for an arctangent test function and the uniform distribution at different maximum polynomial degrees pmax are given in Figure 3 for one random parameter nξ = 1. The error at a constant polynomial degree p = pmax shows a linear convergence of approximately the order O = (p + 1)/nξ = p + 1 for nξ = 1. Until the number of samples ns is sufficiently large to construct an approximation of the maximum degree pmax , the highest possible degree p is used following (15) and (16). This results for pmax = 9 in a continuously increasing polynomial degree and a spectral character of the convergence until it reaches machine precision after ns = 10 samples. Limiting the polynomial degree p only by the number of samples ns and not by a maximum degree pmax results, therefore, in the lowest error and a superlinear convergence. Figure 3b shows that error estimate εˆL2 (25) gives an accurate conservative estimate of both the linear convergence for pmax = 1 and the superlinear convergence for pmax = 9. 0
0
10
10 10
10
−4
10
−4
10
p=1 p=2
−6
p=3 −8
p=4 p=5 p=6 p=7 p=8
10
−10
10
p=1
10
error
−6
10
error
SESC estimate
−2
−2
−8
10
−10
10
−12
−12
10
10
p=9
p=9
−14
10
−14
0
10
1
10
2
10 samples
10
0
1
10
(a) Error for p = {1, . . . , 9}
10 samples
2
10
(b) Error estimate for p = {1, 9}
Figure 3. Randomized refinement sampling for nξ = 1 and the uniform distribution.
Randomized refinement sampling is compared to using higher degree Newton–Cotes quadrature22 in Figure 4 up to nξ = 3 random parameters. Owing to the superlinear convergence, randomized refined sampling converges faster than the linear second–order behavior of Newton–Cotes quadrature. This leads to a reduction of the number of samples ns up to two orders of magnitude for obtaining a same accuracy. Error estimate (25) accurately predicts the error and its dependence on nξ in Figure 5 for nξ = {1, . . . , 4}. The initial error is overpredicted, since for the samples of the initial grid the error is taken to be equal to εk = −vk for k = 1, . . . , nsinit . The asymptotic average number of samples per element ns /ne decreases with increasing dimensionality nξ to a value significantly lower than unity. Table 1 gives the decreasing value of ns /ne up to ns /ne = 0.04 for nξ = 4 at the maximum number of samples used in Figure 5. The number of samples in the initial discretization nsinit is also significantly lower than for Newton–Cotes sampling in Table 1. However, it still increases exponentially with dimensionality nξ following nsinit = 2nξ + 1. This exponential increase of nsinit with nξ is caused by sampling the vertices of the hypercube probability space Ξ. Table 1. Asymptotic average number of samples per element ns /ne and initial number of samples nsinit for randomized refinement (RR) sampling and Newton–Cotes (NC) quadrature. ns ne
nξ 1 2 3 4
NC 1.00 4.59 22.64 -
RR 1.14 0.51 0.16 0.04
nsinit NC RR 3 3 25 5 341 9 17
7 of 17 American Institute of Aeronautics and Astronautics
0
10
SESC RR SESC NC
−2
10
−4
error
10
−6
10
−8
1
10
2 1
3
2
3
−10
10
0
10
1
2
10
10
3
samples
4
10
5
10
10
Figure 4. Randomized refinement (RR) sampling and higher degree Newton–Cotes (NC) quadrature for nξ = {1, 2, 3} and the uniform distribution.
0
10
SESC estimate
−2
10
−4
error
10
−6
10
4
−8
10
1
2
3
−10
10
0
10
1
10
2
samples
10
3
10
Figure 5. Error estimate for randomized refinement sampling with nξ = {1, . . . , 4} and the uniform distribution.
8 of 17 American Institute of Aeronautics and Astronautics
IV.
Essentially Extremum Diminishing extrapolation
The SESC method is combined with extrapolation to achieve a linear increase of the initial number of samples nsinit with nξ . This extrapolation changes the robustness of the method to Essentially Extremum Diminishing23 (EED) in probability space. Refinement criteria for this formulation are applied to the beta and normal distribution. The flexibility of both simplex elements and extrapolation for discretizing complex geometries is also employed to handle non–hypercube probability spaces. Finally the resulting SESC algorithm is summarized. A.
SESC with extrapolation
In order to avoid the exponential increase of the initial number of samples nsinit with nξ , the vertices of the hypercube probability space Ξ are not used as samples in the initial discretization. This introduces the need for extrapolation to approximate the response surface u(ξ) in the entire parameter space Ξ. The vertices of Ξ are still used in constructing the Delaunay triangulation, see Figure 6a for a two–dimensional example. However, the samples vk in the vertices ξ k of Ξ with k = 1, . . . , 2nξ are not computed as denoted by the open circles. If an element Ξj has at least one vertex ξ kj,l for which holds kj,l ≤ 2nξ with l = 0, . . . , nξ , then approximation wj (ξ) in Ξj is constructed using extrapolation. These extrapolation elements are depicted by the dotted element boundaries in Figure 6a. The approximation w(ξ) in the extrapolation elements is found by extending the interpolations wj (ξ) from the other elements Ξj to the boundary ∂Ξ of Ξ. The function w(ξ MCk ) used in a Monte Carlo integration point ξ MCk in an extrapolation element is then equal to wj (ξ) in the nearest interpolation element Ξj .
0.8
0.8
0.6
0.6
ξ2
1
ξ2
1
0.4
0.4
0.2
0.2
0 0
0.2
0.4
ξ
0.6
0.8
0 0
1
0.2
1
0.4
ξ
0.6
0.8
1
1
(a) Triangulation
(b) Extrapolation
Figure 6. SESC discretization of two–dimensional parameter space Ξ with nξ = 2 for EED extrapolation.
S Consider extrapolation domain Ξex = j Ξj for all j = 1, . . . , ne for which kj,l ≤ 2nξ for at least one l ∈ {0, . . . , nξ }. Interpolation region Ξin is then given by Ξ = Ξin ∪ Ξex and Ξin ∩ Ξex = ∅. Consider also the Monte Carlo integration points ξ MCk ∈ Ξex in the extrapolation domain for k = 1, . . . , nsMCex with nsMCex < nsMC and ξMCk ∈ Ξin for k = nsMCex + 1, . . . , nsMC in the interpolation elements. Then the approximation w(ξ MCk ) in the extrapolation domain ξ MCk ∈ Ξex is given by w(ξ MCk ) = wj (ξ MCk ),
(26)
with j minimizing the Euclidean distance min kξMCk − ξ MCl k l
for ξMCl ∈ Ξj ⊂ Ξin ,
(27)
with k = 1, . . . , nsMCex and l = nsMCex + 1, . . . , nsMC . This leads to the extension of the functions wj (ξ) from the interpolation elements Ξj ∈ Ξin to the boundary of Ξ along to the dotted lines in Figure 6b. These regions of extrapolation do not coincide with the boundary elements Ξj ⊂ Ξex of Figure 6a. They form a more suitable basis for extrapolation than the skewed boundary elements of Figure 6a.
9 of 17 American Institute of Aeronautics and Astronautics
The initial discretization contains in this case one sample in the interior of Ξ. The first approximation is obtained after refinement of the grid to nξ + 1 samples, for which the initial linear interpolation can be constructed. This number of samples increases linearly with dimensionality. The element refinement measure that is used for randomized refinement sampling and extrapolation is ¯ 2Oj ¯ j Ξj ε¯L2 j = Ω , (28) ¯ Ξ ¯ needs to be included, since Oj can vary from element to element according to where Ξ Oj =
pj + 1 , nξ
(29)
with pj in element Ξj possibly limited by the LEC condition (17). Refinement measure (28) has a clear relation with moment equation (11) ne Z X µui (x, t) ≈ µwi (x, t) = wj (x, t, ξ)i fξ (ξ)dξ. j=1
Ξj
¯ j reflects the weighting of the integral by density fξ (ξ). The contribution of the error in the Probability Ω ¯ j /Ξ) ¯ 2Oj . The factor two originates from the approximation of wj (x, t, ξ)i is approximated by the term (Ξ L2 error norm. The dependence of local relative error estimate (28) on pj makes ε¯L2 j an indirect solution– based refinement criterion derived from error estimate εˆL2 (25). It is also a weighted refinement measure ¯ j . This weighting becomes less important for increasing order of for non–uniform distributions through Ω convergence Oj . In case of multiple outputs of interest u(ξ), refinement measure ε¯L2 j can be based on either the average of the values of pj or polynomial degree p from (15) and (16). The degree pj for (29) in the extrapolation elements Ξj ⊂ Ξex is determined by averaging the polynomial degree of the approximations w(ξ MCk ) for all ξ MCk ∈ Ξj . Error estimate εˆj in the extrapolation elements Ξj ⊂ Ξex is also computed as the average over Ξj of the estimates εˆl (22) for the approximations w(ξ MCk ) = wl (ξ MCk ) with Ξl ⊂ Ξin and ξ MCk ∈ Ξj . B.
Linear increase of initial number of samples
Results of SESC with extrapolation are compared to the ED method of section III in terms of accuracy, initial number of samples, and robustness. Results of EED extrapolation and the ED formulation of SESC are shown in Figure 7 for the uniform distribution and nξ = {1, . . . , 4}. The EED formulation leads initially to a faster convergence and a lower error for multiple random parameters, e.g. nξ = {3, 4}, owing to the better spread of the samples. Asymptotically, both methods give a comparable error, which suggests that the extrapolation has no significant negative effect on the accuracy. The linear increase of the initial number of samples nsinit of EED extrapolation is significantly slower than the exponential increase of the ED formulation given in Table 2 up to nξ = 8 random parameters. The initial number of samples of the ED method reaches nsinit = 257 for nξ = 8 compared to nsinit = 9 for EED extrapolation. The errors for the initial discretizations of the EED formulation are also given in Table 2. The L2 error norm εL2 is consistently larger than the errors for the mean εµ and standard deviation εσ , which confirms that it is a suitable measure to be used as refinement stopping criterion. The errors of Table 2 further converge by refining the initial discretization to a larger number of elements. ¯ EED in terms of the The robustness of EED extrapolation is quantified in Figure 8a using measure Ω upper bound of the probability of unphysical predictions for nξ = {1, 2, 3}. The ED formulation of SESC ¯ EED = 0. EED extrapolation converges to this value up results by definition in the optimal value of Ω −1 ¯ to ΩEED = 10 for the considered sample sizes. The robustness of the EED formulation can further be ¯ EED . This can, for example, be achieved by increasing the improved by accelerating the convergence of Ω refinement measure ε¯L2 j (28) in the extrapolation elements Ξj ⊂ Ξex by a factor cEED ≥ 1 ε¯L2j
¯j = cEED Ω
¯ 2Oj Ξj , ¯ Ξ
(30)
¯ EED with increasing cEED if Ξj ⊂ Ξex . Results for cEED = {1, 2, 5, 10, 1000} show a faster convergence of Ω ¯ EED = 10−2 in Figure 8b. The convergence also becomes smoother for increasing cEED and nξ , up to Ω 10 of 17 American Institute of Aeronautics and Astronautics
0
10
SESC EED SESC ED
−2
10
−4
error
10
4
−6
10
−8
10
1 3
2
−10
10
0
1
10
10
2
samples
3
10
10
Figure 7. Error convergence for SESC with (EED) and without (ED) extrapolation for nξ = {1, . . . , 4} and the uniform distribution.
Table 2. Initial number of samples nsinit for SESC with (EED) and without (ED) extrapolation, and the corresponding EED errors.
ED nξ 1 2 3 4 5 6 7 8
nsinit 3 5 9 17 33 65 129 257
EED nsinit 2 3 4 5 6 7 8 9
εµ 3.047 · 10−3 2.269 · 10−2 2.346 · 10−2 2.195 · 10−2 2.030 · 10−2 2.263 · 10−2 2.495 · 10−2 2.565 · 10−4
εσ 2.903 · 10−3 2.229 · 10−2 1.058 · 10−1 5.274 · 10−1 3.267 · 10−2 1.160 · 10−1 4.982 · 10−2 2.392 · 10−1
εL2 5.033 · 10−3 4.198 · 10−2 3.383 · 10−1 6.497 · 10−1 5.405 · 10−2 1.672 · 10−1 9.429 · 10−2 3.204 · 10−1
11 of 17 American Institute of Aeronautics and Astronautics
because it concentrates refinement in Ξex . The value of cEED does not significantly affect the approximation accuracy εL2 . 0
0
10
10
1, 2
EED
EED
1 −1
10
2
−1
10
1
5, 10, 1000
3
1
2
2
5
5 10 10 −2
10
1000
−2
0
1
10
10
2
samples
10
3
10
10
0
1
10
10
¯ EED (a) Measure Ω
2
samples
1000 3
10
10
(b) cEED = {1, 2, 5, 10, 1000}
¯ EED for EED extrapolation with the uniform distribution and nξ = {1, 2, 3}. Figure 8. Robustness measure Ω
C.
Non–uniform distributions
Refinement measure ε¯L2 j (28) is in this section applied to the beta and normal distribution. The unbounded parameter domain of the normal distribution is treated by truncating Ξ beyond the Monte Carlo integration point ξMCk that is located furthest in the tails of the normal distribution, which does not affect the approximation of the statistical moments. The error for the beta and normal distribution is compared to the results for the uniform distribution in Figure 9a for cEED = 1. The beta distribution results in a lower error owing to the lower input standard deviation than the uniform distribution in the same domain Ξ. The larger domain Ξ for the normal distribution with the same standard deviation as the uniform distribution results in a slower error convergence. The accuracy depends, therefore, more on the input standard deviation and the size of domain Ξ, than on ¯ EED also decreases signifwhether the input distribution is uniform or non–uniform. Robustness measure Ω icantly faster for the beta and normal distribution in Figure 9b up to an upper bound for the probability ¯ EED = 2.1 · 10−5 for the beta distribution with nξ = 3. This is caused by the of unphysical predictions of Ω combination of the decreasing probability density in the tails of the distribution near the boundary of Ξ and refinement measure ε¯L2 j . 0
0
10
10
−2
10
1
−1
10
2
3
−4
10
1 1
−2
−6
EED
error
10 3
10
1
−3
10
2
2
−8
10
1
10
uniform beta normal
−12
10
0
10
−4
10
2
1
2 1
10
2
samples
10
3
2
3 −10
uniform beta normal
3 3
0
10
10
3 1
10
2
samples
10
3
10
¯ EED (b) Robustness measure Ω
(a) Error
¯ EED for EED extrapolation with the uniform, beta, and normal Figure 9. Error and robustness measure Ω distribution for nξ = {1, 2, 3}.
12 of 17 American Institute of Aeronautics and Astronautics
The resulting discretizations of Ξ for the beta and normal distribution are given in Figure 10. It shows the effect of refinement measure ε¯L2 j in terms of a good coverage of domain Ξ by the elements combined with a moderately higher density of samples at higher probabilities in the center of the domain. 1 1.5 0.8 1
ξ2
ξ2
0.6 0.5
0.4 0 0.2 −0.5 0 0
0.2
0.4
ξ
0.6
0.8
1
−0.5
0
0.5 ξ
1
1.5
1
1
(a) Beta distribution
(b) Normal distribution
Figure 10. EED extrapolation discretizations for the beta and normal distribution in a two–dimensional parameter space Ξ with nξ = 2 and ns = 136.
D.
Non–hypercube probability spaces
The problems considered so far have in common that the parameter domain Ξ is assumed to be a hypercube. However, in practice the parameter domain can have other shapes if the random input parameters are correlated. For certain regions of the parameter space either the probability can be zero or the physical system can even be unrealizable. This can occur, for example, in case of uncertainty originating from production tolerances, where products with a too high deviation from the design geometry are rejected. This leads to truncation of the probability space, which can be irregular for multi–dimensional parameterizations of the geometrical variability. The flexibility of simplex elements with extrapolation is used here to discretize these non–hypercube probability spaces. Consider a non–hypercube parameter domain Ξ illustrated in Figure 11 by the gray area. Assume that the deterministic problem does not have a physical solution for the parameter values ξ outside domain Ξ. The samples vk for the sampling points outside the parameter domain ξk ∈ / Ξ can in that case not be computed as denoted by the open circle. These sampling points are treated similarly as the vertices of the / Ξ are used in constructing the hypercube probability space in section IV.A. The sampling points ξ k ∈ Delaunay triangulation, however, they are not used for building the approximation w(ξ) since the samples vk are not available. Approximation wj (ξ) in elements with at least one of its vertices ξ kj,l with l = 0, . . . , nξ outside Ξ denoted by the dotted lines in Figure 11, is determined using extrapolation in the same way as for the extrapolation elements Ξj ⊂ Ξex as described in section IV.A. Elements Ξj ⊂ Ξin with all vertices ξ kj,l inside Ξ determine wj (ξ) using interpolation and a higher degree stencil Sj as before, where Ξin does in general not coincide with Ξ. The initial discretization consists of points in the vertices of the hypercube enclosing Ξ and one sample in its interior. The error estimates and the refinement procedure of EED extrapolation of section IV.A are ¯ j and probability Ω ¯ j of the elements Ξj the region inside Ξ is considered only. used, where for the volume Ξ An example of a non–hypercube probability space Ξ is shown in Figure 12 with its MC and SESC discretizations. The probability inside Ξ is uniformly distributed as illustrated by the MC sampling of Figure 12a. In this case the random parameters ξ = (ξ1 , ξ2 ) are correlated through the boundary ∂Ξ of Ξ. Refinement criterion ε¯L2 j results in an effective simplex elements discretization of the interior of domain Ξ in Figure 12b. The points ξ k outside Ξ and the extrapolation elements Ξj ⊂ Ξex are denoted by open circles and dotted lines, respectively. The dots and the bold lines show the samples and the interpolation elements Ξj ⊂ Ξin inside Ξ. A second example is given in Figure 13 where the lower left corner is excluded from square domain Ξ. This example shows that it is possible to combine sampling of the vertices of the square Ξ at (ξ1 , ξ2 ) equal to (0, 1), (1, 0), and (1, 1), with extrapolation to the truncation boundary of Ξ in the 13 of 17 American Institute of Aeronautics and Astronautics
Figure 11. Example discretization of a two–dimensional non–hypercube probability space Ξ denoted by the gray area.
lower left corner. The method can also be used for non–uniform probability densities in the non–hypercube domains Ξ and for more than two random parameters.
0.8
0.8
0.6
0.6
ξ2
1
ξ2
1
0.4
0.4
0.2
0.2
0 0
0.2
0.4
ξ
0.6
0.8
0 0
1
0.2
1
0.4
ξ
0.6
0.8
1
1
(a) MC with nsMC = 104
(b) SESC with ns = 40
Figure 12. MC and SESC discretization of a two–dimensional non–hypercube probability space Ξ with nξ = 2.
The error convergence for the two non–hypercube probability domains Ξ of Figures 12 and 13 is given in Figure 14. The non–hypercube shape of domain Ξ does not affect the superlinear character of the convergence. The error is even slightly smaller than for the hypercube Ξ and the uniform distribution of Figure 7, owing to the smaller size of Ξ in the non–hypercube examples. E.
SESC algorithm summary
The SESC algorithm including randomized refinement sampling and EED extrapolation can be summarized as follows for hypercube and non-hypercube probability spaces: 1. The initial grid of sampling points ξ k is composed of the 2nξ vertices of the hypercube enclosing the probability space Ξ and one sampling point in the interior. 2. The nsinit samples vk are computed by solving deterministic problems (4) only for the nsinit parameter values of the initial sampling points ξk located in Ξ. 3. Polynomial order p is determined from the available number of samples ns using (15) and (16).
14 of 17 American Institute of Aeronautics and Astronautics
0.8
0.8
0.6
0.6
ξ2
1
ξ2
1
0.4
0.4
0.2
0.2
0 0
0.2
0.4
ξ
0.6
0.8
0 0
1
0.2
0.4
ξ
0.6
0.8
1
1
1
(a) MC with nsMC = 104
(b) SESC with ns = 40
Figure 13. MC and SESC discretization of a two–dimensional non–hypercube probability space Ξ with nξ = 2.
0
10
domain 1 domain 2
−2
10
−4
error
10
−6
10
−8
10
−10
10
0
10
1
10
2
samples
10
3
10
Figure 14. Error convergence for the non–hypercube probability domains 1 and 2 of Figures 12 and 13, respectively.
15 of 17 American Institute of Aeronautics and Astronautics
4. The Delaunay triangulation consisting of ne simplex elements Ξj is constructed based on all sampling points ξ k . 5. The extrapolation elements Ξj ⊂ Ξex are identified for which hold for at least one of the vertices ξkj,l ∈ / Ξ with l = 0, . . . , nξ . The other elements form the interpolation region Ξin , with Ξ = Ξin ∪ Ξex . 6. The polynomial approximation wj (ξ) (8) in each of the interpolation elements Ξj ⊂ Ξin is constructed in the following steps: a) The stencil Sj of Nj + 1 samples vj for polynomial order pj = p is selected based on (13) with Nj from (18). b) The polynomial coefficients cj,l of wj (ξ) for l = 0, . . . , Nj are computed by solving (10). c) If interpolation wj (ξ) is not LEC in element Ξj with respect to the samples vj by violating (17), then local polynomial degree pj is decreased by one and the loop for element Ξj is returned to step a) to select a reduced stencil Sj . Interpolation wj (ξ) for pj = 1 is LEC by definition. d) Local error estimate εˆj (22) and refinement measure ε¯L2 j (28) are computed. 7. Approximation wj (ξ) in the other elements Ξj ⊂ Ξex is found using extrapolation according to (26) and (27). Local error estimate εˆj and refinement measure ε¯L2 j in the extrapolation elements Ξj ⊂ Ξex are determined from averaging over Ξj as detailed in section IV.A. 8. The probability distribution and the statistical moments (11) of w(ξ) are computed using Monte Carlo integration, for example, for the mean µw and standard deviation σw . The local error estimates εˆj are also combined into global error estimate εˆL2 (25). 9. If global error estimate εˆL2 is smaller than a user–defined threshold value the algorithm is stopped. 10. Otherwise, the element Ξj with the largest value of refinement measure ε¯L2 j is refined by randomly selecting a new sampling point ξk in sub–simplex Ξsubj defined by (12). The algorithm can be parallelized by refining multiple elements neref ≥ 1 in one refinement step. 11. If the new sampling point ξ k is located in Ξ, then deterministic problem (4) is solved for the parameter values ξk to compute the corresponding sample value vk . 12. The algorithm finally returns to step 3 to construct a new approximation w(ξ) for the new set of sampling points ξ k . In case of a hypercube probability space Ξ, the ED formulation of SESC from section III is retrieved by considering Ξ including its boundary ∂Ξ in steps 2 and 11. Considering Ξ without its boundary ∂Ξ leads to SESC with EED extrapolation from section IV.A.
V.
Conclusions
The Simplex Elements Stochastic Collocation (SESC) method is introduced for robust and efficient propagation of multiple uncertainties in computational models. The non–intrusive Uncertainty Quantification (UQ) method is equally robust as Monte Carlo (MC) simulation in terms of the Extremum Diminishing (ED) concept in probability space. SESC combines robustness with the efficiency of the superlinear convergence of Stochastic Collocation (SC) methods and a linear increase of the initial number of samples with dimension. The efficiency of the simplex elements discretization of probability space is based on higher degree polynomial interpolation, randomized refinement sampling, and Essentially Extremum Diminishing (EED) extrapolation. Randomly selecting the new sampling point at the refinement of a simplex leads to superlinear convergence and a reduction of samples up to two orders of magnitude compared to using higher–degree Newton–Cotes quadrature points in the numerical results for an analytical test function verification. The combination of SESC with extrapolation results in the linear increase of the initial number of samples with nξ . The robustness is in that case quantified by the upper bound for the probability of unphysical predictions ¯ EED by introducing the EED concept. Ω The adaptive refinement of SESC employs a solution–based refinement measure ε¯L2 j based on a local relative error estimate to refine one element Ξj at a time or more in parallel. All samples vk are reused after 16 of 17 American Institute of Aeronautics and Astronautics
the nested refinements, until a global absolute error estimate εˆL2 reaches a refinement stopping criterion. The asymptotic average number of samples per element ns /ne decreases with increasing nξ to a value significantly lower than unity of, for example, ns /ne = 0.04 for nξ = 4. The application of the flexibility of SESC to examples of non–hypercube probability spaces Ξ with correlated random parameters results in a faster convergence than for their enclosing hypercube probability space. The focus of the current paper is the efficient extension of the simplex elements discretization to multiple random parameters. In future work the effectiveness of the multi–element nature of SESC for solving discontinuous and non–monotonic problems will be examined in more detail for computational applications.
Acknowledgements This work was supported by the United States Department of Energy under the Predictive Science Academic Alliance Program (PSAAP) at Stanford University.
References 1 I. Babuˇ ska, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal. 45 (2007) 1005–1034. 2 M. De Berg, O. Cheong, M. Van Kreveld, M. Overmars, Computational Geometry: Algorithms and Applications, Springer–Verlag, Berlin (2008). 3 M. Deb, I. Babuˇ ska, J. Oden, Solution of stochastic partial differential equations using Galerkin finite element techniques, Comput. Method. Appl. M. 190 (2001) 6359–6372. 4 J. Foo, X. Wan, G.E. Karniadakis, The multi–element probabilistic collocation method (ME-PCM): Error analysis and applications, J. Comput. Phys. 227 (2008) 9572–9595. 5 J. Foo, G.E. Karniadakis, Multi–element probabilistic collocation method in high dimensions, J. Comput. Phys. 229 (2010) 1536–1557. 6 B. Ganapathysubramanian, N. Zabaras, Sparse grid collocation schemes for stochastic natural convection problems, J. Comput. Phys. 225 (2007) 652–685. 7 R.G. Ghanem, P. Spanos, Stochastic finite elements: a spectral approach, Springer–Verlag, New York (1991). 8 J.M. Hammersley, D.C. Handscomb, Monte Carlo Methods, Fletcher & Son Ltd, Norwich (1964). 9 S. Hosder, R.W. Walters, R. Perez, A non–intrusive polynomial chaos method for uncertainty propagation in CFD simulations, 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada (2006) AIAA-2006-891. 10 A. Jameson, Positive schemes and shock modelling for compressible flows, Int. J. Num. Meth. Fluids 20 (1995) 743–776. 11 X. Ma, N. Zabaras, An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations, J. Comput. Phys. 228 (2009) 3084–3113. 12 O.P. Le Maˆ ıtre, O.M. Knio, H.N. Najm, R.G. Ghanem, Uncertainty propagation using Wiener-Haar expansions, J. Comput. Phys. 197 (2004) 28–57. 13 O.P. Le Maˆ ıtre, H.N. Najm, R.G. Ghanem, O.M. Knio, Multi–resolution analysis of Wiener–type uncertainty propagation schemes, J. Comput. Phys. 197 (2004) 502–531. 14 H.N. Najm, Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics, Annu. Rev. Fluid Mech. 41 (2009) 35–52. 15 F. Nobile, R. Tempone, C.G. Webster, A sparse grid stochastic collocation method for partial differential equations with random input data, SIAM J. Numer. Anal. 46 (2008) 2309–2345. 16 M.A. Tatang, Direct incorporation of uncertainty in chemical and environmental engineering systems, PhD thesis, MIT, Cambridge (1995). 17 X. Wan, G.E. Karniadakis, Multi–element generalized polynomial chaos for arbitrary probability measures, SIAM J. Sci. Comput. 28 (2006) 901–928. 18 N. Wiener, Nonlinear problems in random theory, MIT Technology Press and John Wiley & Sons, New York (1958). 19 J.A.S. Witteveen, G.J.A. Loeven, H. Bijl, An adaptive stochastic finite elements approach based on Newton–Cotes quadrature in simplex elements, Comput. Fluids 38 (2009) 1270–1288. 20 J.A.S. Witteveen, H. Bijl, Effect of randomness on multi–frequency aeroelastic responses resolved by unsteady adaptive stochastic finite elements, J. Comput. Phys. 228 (2009) 7025–7045. 21 J.A.S. Witteveen, H. Bijl, A TVD uncertainty quantification method with bounded error applied to transonic airfoil flutter, Commun. Comput. Phys. 6 (2009) 406–432. 22 J.A.S. Witteveen, G. Iaccarino, Simplex elements stochastic collocation for uncertainty propagation in robust design optimization, 48th AIAA Aerospace Sciences Meeting, Orlando, Florida (2010) AIAA-2010-1313. 23 J.A.S. Witteveen, G. Iaccarino, Simplex elements stochastic collocation (2010) submitted. 24 D. Xiu, G.E. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput. 24 (2002) 619–644. 25 D. Xiu, J.S. Hesthaven, High–order collocation methods for differential equations with random inputs, SIAM J. Sci. Comput. 27 (2005) 1118–1139. 26 D. Xiu, Fast numerical methods for stochastic computations: A review, Commun. Comput. Phys. 5 (2009) 242–272.
17 of 17 American Institute of Aeronautics and Astronautics
