Minimal multi-element stochastic collocation for uncertainty ...

Journal of Computational Physics 242 (2013) 790–808

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions John D. Jakeman a,⇑,1, Akil Narayan b,2, Dongbin Xiu c,2 a

Sandia National Laboratories, Albuquerque, NM 87185, United States Department of Mathematics, University of Massachusetts Dartmouth, 285 Old Westport Road, North Dartmouth, MA 02747, United States c Department of Mathematics, Purdue University, West Lafayette, IN 47907, United States b

a r t i c l e

i n f o

Article history: Received 22 February 2012 Received in revised form 12 February 2013 Accepted 13 February 2013 Available online 13 March 2013 Keywords: Uncertainty quantification Generalized polynomial chaos Stochastic collocation Multi-element Discontinuous functions

a b s t r a c t We propose a multi-element stochastic collocation method that can be applied in highdimensional parameter space for functions with discontinuities lying along manifolds of general geometries. The key feature of the method is that the parameter space is decomposed into multiple elements defined by the discontinuities and thus only the minimal number of elements are utilized. On each of the resulting elements the function is smooth and can be approximated using high-order methods with fast convergence properties. The decomposition strategy is in direct contrast to the traditional multi-element approaches which define the sub-domains by repeated splitting of the axes in the parameter space. Such methods are more prone to the curse-of-dimensionality because of the fast growth of the number of elements caused by the axis based splitting. The present method is a two-step approach. Firstly a discontinuity detector is used to partition parameter space into disjoint elements in each of which the function is smooth. The detector uses an efficient combination of the high-order polynomial annihilation technique along with adaptive sparse grids, and this allows resolution of general discontinuities with a smaller number of points when the discontinuity manifold is low-dimensional. After partitioning, an adaptive technique based on the least orthogonal interpolant is used to construct a generalized Polynomial Chaos surrogate on each element. The adaptive technique reuses all information from the partitioning and is variance-suppressing. We present numerous numerical examples that illustrate the accuracy, efficiency, and generality of the method. When compared against standard locally-adaptive sparse grid methods, the present method uses many fewer number of collocation samples and is more accurate. Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction Uncertainty quantification (UQ) and stochastic computation have received much attention in recent years, with many research efforts devoted to the development of efficient numerical methods for complex systems. A particular focus is on the design of algorithms that are more efficient than the traditional Monte Carlo (MC) method, which, though extremely robust and easy to implement, can be prohibitively time-consuming for practical systems when high accuracy is required.

⇑ Corresponding author. Tel.: +1 5052849097. E-mail addresses: [email protected] (J.D. Jakeman), [email protected] (A. Narayan), [email protected] (D. Xiu). Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the US Department of Energys National Nuclear Security Administration under contract DE-AC04-94AL85000. 2 This was supported by AFOSR, DOE/NNSA, and NSF. 1

0021-9991/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.02.035

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

791

One of the more widely adopted methodologies is generalized polynomial chaos (gPC) [33], an extension of the standard polynomial chaos (PC) method [17]. Here one seeks to construct accurate polynomial approximations of the unknown target function in the underlying random space. Whenever the target function is sufficiently smooth, the gPC approximation becomes highly accurate even at low orders and the methodology is highly efficient. To construct such gPC polynomial approximations, the traditional approach is to employ stochastic Galerkin (SG), where the residual of the stochastic equations is orthogonalized against the polynomial space [6,17,31,33]. Though mathematically rigorous and accurate, SG method typically results in a large system of coupled deterministic equations and requires one to develop new simulation codes. An alternative approach, widely used in many practical systems after its high-order formulation [32], is stochastic collocation (SC). Here one conducts sampling-like simulations, just like MC, and then constructs polynomial approximations at the post-processing step [5,24,28,32]. Therefore, SC combines the strengths of non-intrusive sampling methods such as MC based methods and the fast convergence of gPC-SG methods. The performance of both SG and SC methods depend critically on the smoothness of the target function. Whenever the function is smooth in the random parameter space, these methods converge quickly and are highly efficient. On the other hand, their performance deteriorates when the target function lacks regularity, and in particular, possesses discontinuities. In these situations, the gPC methods using global polynomials suffer from Gibbs-type phenomenon and have very slow convergence. To circumvent the difficulty, multi-element gPC (ME-gPC) methods were developed [29,30]. The idea is to decompose the random space into sub-domains, in each of which the target function is smooth and amenable to local gPC approximations using either SG or SC. The final global approximation in the entire domain is then a ‘‘patched’’ solution of local solutions in all elements [1,8,12,13,19,22,23]. The undesirable impact of the discontinuities is thus confined in a limited number of elements surrounding the discontinuities, and the global solution can regain high accuracy away from them. The challenge of the standard multi-element (ME) approaches is its simulation cost. One now needs to first resolve the stochastic problem in each of the elements. In the existing ME methods, the elements are constructed by splitting each axis, and then defining the corresponding hypercubes. A prominent drawback of this approach is that the construction inevitably utilizes a tensor structure, which results in a fast growth of the total number of elements in high dimensional random space. For example, if each axis is split into two parts (the minimal amount of splitting), then the total number of elements in ddimensional random space is 2d , where each element requires a solution of the original d-dimensional stochastic system. In high dimensions d  1, the total simulation cost can be prohibitive. (In many cases, the axes are required to be split into more than two parts, though adaptive algorithms can reduce the number of splitting.) In this paper we present a different multi-element method that abandons the tensor structure in local element constructions. A distinct feature of the current method is that the local elements are defined by the underlying stochastic problem directly. More precisely, we seek to decompose the random space by splitting it into elements along where the discontinuities lie. By doing so, the total number of elements equals the number of smooth sub-domains defined by the underlying target function. In this sense the current method can be considered optimal. Hence the term minimal-element (mE) method. The minimal-element method consists of the following important components. First, we need an algorithm to efficiently detect the existence of discontinuities; and if there is one, to locate its geometry and classify the sub-domains containing smooth parts of the function. In the sequel, we will refer to sub-domains that contain smooth parts of the function as smooth sub-domains. To this end, we employ the high-order polynomial edge detection method using adaptive sparse grids [20]. This method is based on the polynomial annihilation edge detection method in one and two dimensions [3], which was shown to hold certain advantages over some of other more traditional methods ([10,14,27]). It was first generalized to high dimensional spaces via a straightforward dimension-by-dimension extension [2]. The work of [20] significantly improves the performance of the method by using an adaptive sparse grids algorithm. Once the random space has been decomposed into disjoint elements, defined by the smooth sub-domains of the target function, the next task is to construct accurate polynomial approximations in each elements. The challenge here is that the elements are of irregular shape, because of the arbitrary geometry of the underlying discontinuities. Attempts to map the irregular elements into regular elements will be difficult, if not impossible, especially in high dimensional spaces. Our proposed strategy is to use the stochastic collocation (SC) method directly on the irregular elements. In particular, since we have already computed the solution ensembles in the discontinuity detection step, we will not conduct further SC simulations. Instead, we will seek to construct gPC approximations in each element using the existing simulation results on the sparse grid generated by the discontinuity detector. The difficulty is that now the collocation points do not possess any structure because of the arbitrary boundaries imposed by the discontinuities. Also the number of collocation points in each element can be arbitrary. In order to construct high-order gPC polynomial approximations in each element, we employ the ‘‘least orthogonal interpolation’’ method developed in [25]. This method allows one to construct high order polynomial approximations in high dimensions based on arbitrary number of collocation nodes located at arbitrary locations. In this paper, we employ this method and conduct a further improvement. We adaptively choose subsets of the sparse grids, from the discontinuity detection step, and construct an interpolant that minimizes oscillations. The result is a nonlinear interpolation method, robust for smooth functions, that can perform post-processing polynomial construction regardless of nodal distribution or Euclidean dimension. We show that this method performs well when applied to point sets given by the discontinuity detector. We remark that the least orthogonal interpolation is merely a choice we make here. One is free to use other technique for the polynomial approximation in each element. For example, one can employ a least-square type polynomial regression. The least orthogonal interpolation is employed here largely due to many of its desirable mathematical properties

792

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

[25]: its formation is unique given a nodal distrubtion and a weight function; and the interpolant is continuous, within reason, with respect to the weight function and the nodal distribution. This paper is outlined as follows: Section 2 gives the general setup for a multi-element stochastic method. Section 3 reviews some of the theory behind the polynomial annihilation detection technique and its application in the context of highdimensional adaptive sparse grids. Section 4 provides a brief summary of the least orthogonal interpolant and introduces the adaptive scheme we use for reconstruction. Finally, Section 5 produces examples that show the efficiency and accuracy of our method when applied to test problems. 2. Multi-element methods for stochastic problems Our goal is to solve a general stochastic problem

Lðf ; xÞ ¼ 0; where x ¼ ðx1 ; . . . ; xd Þ 2 Rd ; d P 1, is the multidimensional random parameter characterizing the uncertain inputs, and f is the unknown solution. The operator L can be a nonlinear partial differential equation, in which case we need proper initial and boundary conditions. Here we neglect the explicit description of the spatial and temporal variables and focus only on the solution dependence on the random variables x. In gPC methods, we seek to construct a polynomial

fN ðxÞ  f ðxÞ;

x 2 C # Rd ;

where N typically represents the order of the polynomial approximation and we use C to denote the random space. When f is sufficiently smooth, fN converges to f quickly, and the approximation is highly accurate even at very small order N. However, when f is discontinuous, the convergence of a global approximation deteriorates significantly. It is therefore desirable to use a multi-element (ME) method, where the negative impacts of the discontinuities are confined to a limited number N e of small elements E k surrounding the discontinuities. With this aim we decompose the input space C into N e disjoint elements E k ; k ¼ 1; . . . ; N e such that

C¼

Ne [

Ek;

and E k

\

E j ¼ ;;

if k – j:

ð1Þ

k¼1

Once the random input space has been divided into a number of disjoint elements, a locally defined approximation can be constructed on each element E k from a set of collocation nodes. (Alternatively, stochastic Galerkin can be used in each element as well.) Given a set of approximations IE k ½f ðxÞ on each element E k the global approximation can be constructed by patching the local approximations together,

fNe ðxÞ ¼

Ne X I Ek ½f ðxÞvE k ðxÞ;

ð2Þ

k¼1

where v is the standard indicator function satisfying vA ðxÞ ¼ 1, if x 2 A and zero otherwise. Once again the approximation operator can be constructed via either a Galerkin or a collocation procedure. In term of SC, the usual choices include interpolation, regression, or discrete projection. (For a detailed description of these approximation procedures, see [31].) 2.1. Standard multi-element methods In the existing literature the standard ME methods utilize tensor construction for the element definition. One first decomposes each axis into intervals, where the exact partition depends on the nature of the discontinuities cutting through the axes. For each i ¼ 1; . . . ; d, let ni P 0 be the number of intervals in the xi axis, and ðxi;0 ; . . . ; xi;ni Þ be the coordinates at the interfaces. The elements in the entire space are thus defined as hypercubes based on the tensor products of the one-dimensional coordinates. In particular, the kth-element can be defined as

E k ¼ ðx1;k1 ; x1;k1 þ1 Þ      ðxd;kd ; xd;kd þ1 Þ;

0 6 ki < ni ; i ¼ 1; . . . ; d;

where an ordering scheme is employed so that there is a unique one-to-one correspondence between the single index k and the multi-index ðk1 ; . . . ; kd Þ. Though the general formulation of ME methods applies to any type of element construction, in current practice the hypercube construction using tensorization remains the only choice. A direct consequence of this construction is that the total number of elements is N e ¼ n1      nd , and grows exponentially fast in high dimensions d  1. Considering the fact that in each element one still needs to resolve the original d-dimensional stochastic system, the overall simulation cost can be prohibitively high. Though the size of each element is much smaller than the global domain and adaptive strategies can reduce the total number of element significantly, the tensor construction still remains a quite costly method in high dimensions.

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

793

3. Minimal-element method: domain decomposition In the minimal-element (mE) method, we will not employ the tensor construction for the elements. Instead, we seek to decompose the random space into smooth sub-domains that are defined by the structure of the discontinuities. By doing so the number of elements is minimal because it equals the number of the smooth sub-domains given by the problem. Our domain decomposition method is heavily based on the high-dimensional discontinuity detection algorithm developed in [20]. 3.1. Polynomial annihilation discontinuity detection The discontinuity detection algorithm we utilize employs the one-dimensional polynomial annihilation discontinuity detection method [3] to detect discontinuities on a set of points. The points are found using an adaptive sampling method similar to that used in locally adaptive sparse grid approximation methods [18,21,23]. The polynomial annihilation method is used to guide refinement, thereby investing the majority of function evaluations in regions surrounding discontinuities and neglecting grid points in smooth regions. Consider a piecewise continuous function f ðxÞ : ½0; 1 ! R in one dimensional (random parameter) space that has welldefined one-sided limits, f ðx Þ, at any point x in the domain. Let J denote the set of points of discontinuity of f, that is, J ¼ fx 2 ð0; 1Þ : f ðxþ Þ–f ðx Þg. Let us define the jump function

½f ðxÞ ¼ f ðxþ Þ  f ðx Þ ¼



0;

if x R J; þ



f ðx Þ  f ðx Þ; if x 2 J:

½f  vanishes on ð0; 1Þ n J and is equal to the jump value of f for points in J. Approximating this jump function will allow us to detect and locate any discontinuities in the underlying function f. The polynomial annihilation discontinuity detection method, introduced in [3], seeks an approximation to ½f ðxÞ that converges rapidly to zero away from jump discontinuities. Given a set of unique points Sx ¼ fxj kxj 2 ½0; 1; j ¼ 0; . . . mg at which the value of the function f ðxj Þ is known, the polynomial annihilation method approximates the jump function by

Lm f ðxÞ ¼

1 X cj ðxÞf ðxj Þ; qm ðxÞ x 2S j

ð3Þ

x

The coefficients cj ðxÞ are chosen to annihilate polynomials of degree up to m  1 and are determined by solving the linear system

X

ðmÞ

cj ðxÞpi ðxj Þ ¼ pi ðxÞ;

8 i ¼ 0; . . . ; m:

ð4Þ

xj 2Sx th Here fpi gm derivative of i¼0 is any basis for the space of univariate polynomials of degree m or less, and pi ðxÞ denotes the m pi ðxÞ. The solution to (4) exists and is unique and the value of the coefficients is dependent only on the set Sx and not on the data f ðxj Þ. The normalization factor qm ðxÞ is always non-zero and ensures that Lm f ðxÞ has correct value at the jump discontinuities. The polynomial annihilation discontinuity detection method converges to ½f ðxÞ with a rate depending on m and the local smoothness of f. Specifically, given the maximum separation hðxÞ of the local stencil Sx ðmÞ

hðxÞ ¼ maxfjxi  xi1 j : xi1 ; xi 2 Sx g; and a discontinuity located at f, (3) satisfies the following property [3]:

( Lm f ðxÞ ¼

minðm;kÞ

Oðh

ðxÞÞ;

if f 2 Ck ðIx Þ; k > 0;

½ f ðfÞ þ OðhðxÞÞ; otherwise; where xj1 6 f;

x 6 xj :

ð5Þ

where, Ix is the smallest closed interval such that Sx Ix . The order m in (3) affects the resolution of the jump function both close to and away from any discontinuities. Small m might cause a steep gradient to be misidentified as a discontinuity (due to low resolution). On the other hand, the inherent nature of high order polynomial approximation causes oscillations to occur in the vicinity of the discontinuities when m is large. These oscillations may also be misinterpreted as discontinuities. To prevent inaccuracies due to either of these problems, the minmod function was used in [3] to enhance the performance of the polynomial annihilation method. 3 Specifically, given a set of annhilation orders M ¼ f1; . . . ; #Sx  1g, we compute

8 > < minm2M Lm f ðxÞ; if Lm f ðxÞ > 0; MMðLm f ðxÞÞ ¼ maxm2M Lm f ðxÞ; if Lm f ðxÞ < 0; > : 0; otherwise:

ð6Þ

3 The minmod technique was introduced in the context of edge detection in [7,15], and is typically utilized in flux limiting methods to reduce oscillations when solving numerical conservation laws.

794

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

The minmod function controls the oscillations while still maintaining a high order of convergence away from the jump discontinuities. Examples of its effectiveness can be found in [3]. 3.2. High-dimensional discontinuity detection Extension of the edge detection method to high dimensions requires care because of the potentially exponential growth of the simulation effort. (Most of the existing discontinuity detection methods are restricted to low dimensions d 6 3). Here we use the high-dimensional polynomial annihilation method, proposed in [20], which uses adaptive sparse grids. This method has been shown to enable characterization of discontinuities in function with large numbers of input dimensions. The success of the method arises from its ability to restrict attention to lower-dimensional structure of the discontinuity manifold embedded in the higher-dimensional space. In the case where the function is linear with respect to the input variables not causing the discontinuity, the number of function samples required by the detection method grows linearly with dimension. It must be noted, however, that the number of function evaluations required grows exponentially with the number of dimensions that cause the discontinuity [20]. The adaptive sparse grid is based upon sparse tensor products of nested equidistant one-dimensional grids. Each onedimensional grid point is identified uniquely by a level l and position index i. More specifically

8 l l > < i  2 ; l > 1; 1 < i < 2 ; i odd; l xl;i ¼ i  2 ; l ¼ 1; i 2 f0; 2g; > : 0:5; l ¼ 0: The one-dimensional mesh is shown in Fig. 1. Similarly multi-dimensional sparse grid points nl;i ¼ ðnl1 ;i1 ; . . . ; nld ;id Þ are defined by a multi-dimensional level index l ¼ ðl1 ; . . . ; ld Þ and a multi-dimensional position index i ¼ ði1 ; . . . ; id Þ. The sparse grid is built adaptively. The algorithm begins with an initial low-level sparse grid with all points with P klk1 6 lmin , where typically lmin ¼ 1 and klk1 ¼ di¼1 li . A two-dimensional grid with lmin ¼ 1 is shown in Fig. 2. The value of jump function (3) at the midpoint of each interval (gray points) is then calculated. Any midpoint xl;i with Lm f ðxl;i Þ > e is flagged for refinement (green points). Once all the points that warrant refinement have been determined these points are added to the sparse grid and the children of these points (Gray) are determined. The jump function is then evaluated at each of these children and the refinement process is repeated until a termination condition is met. In one dimension the children of the point xl;i are determined by traversing down the binary tree (refer to Fig. 1). If we denote the parent of a one-dimensional grid point x as PðxÞ, the 2d children of a d-dimensional point x are

fy k ðPðy1 Þ; y2 ; . . . ; yd Þ ¼ x; ðy1 ; Pðy2 Þ; . . . ; yd Þ ¼ x; .. . ðy1 ; y2 ; . . . ; Pðyd ÞÞ ¼ xg where the parent of the root is itself, i.e., Pð0:5; . . . ; 0:5Þ ¼ 0:5. Two steps of the discontinuity detection algorithm are shown in Fig. 2. In the first step the jump function is calculated at all the Gray points by constructing one-dimensional stencils in each dimension. As an example consider the point x2;3;0;0 . To calculate the jump function at this point along the 1st dimension , we must define the local stencil Sx ¼ fx1;0;0;0 ; x0;0;0;0 ; x1;2;0;0 g. Using this stencil we find Lm f ðxl;i Þ > e and thus flag this point for refinement. We cannot calculate the jump function in the second dimension because we need at least two points in the stencil Sx . In the second step the green point, which has been flagged for refinement, is refined. The values of the jump function at the children (Gray) of this point are then calculated. When refining the grid we must ensure that all parents of newly added

Fig. 1. The dyadic splitting of the one-dimension mesh. Children of a point xl;i are the 2 points (in one-dimension, 2d in d-dimensions) with level l þ 1 and connected by one segment in the binary tree.

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

795

Fig. 2. First two steps of the discontinuity detector algorithm. Step 1 (left): Generate initial sparse grid with lmin ¼ 1. (black points), then apply polynomial annihilation to the midpoints (Gray) of all intervals in the initial sparse grid. Step 2 (right): If Lm f ðxl;i Þ > e (green) evaluate the target function at that point and add it to sparse grid. Lastly, apply polynomial annihilation to each of the children of the refined point(s). Repeat step 2 until target resolution d is reached. Note that the function values are known at each (black) point but not the Gray points and for any given midpoint polynomial annihilation can be applied up to d times, once for every dimension. Above each point is labeled according to its level and position multi-indices xðl1 ;l2 Þ;ði1 ;i2 Þ . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

points are already in the sparse grid. For example, when calculating the jump function at x2;3;1;2 we must add any missing parents of that point (red). If any parents are missing they are automatically added to the sparse grid. The discontinuity detection algorithm will continue on indefinitely unless an appropriate termination condition is enforced. With this in mind we a priori set a resolution tolerance d. The choice of d sets the resolution at which any discontinuities present in the target function are resolved. If the distance between a child and its parent is less than d in the ith dimension that point is not refined further in that dimension regardless of the value of the jump function. Due to this restriction, eventually a refinement step will arise during which no points will be flagged for refinement. At this step, the algorithm will terminate. Details of the algorithm, along with pseudocode, are presented in [20]. The method concentrates the vast majority of points in the dimensions that contribute to the discontinuities present in the function. This is achieved by identifying the effective dimension of the manifold and thereafter restricting function evaluations to coordinates that lie on the manifold. Such an approach is unsuited to high-dimensional interpolation. To obtain an accurate interpolant collocation nodes must be generated in regions of variability that do not lie on the discontinuity manifold. To overcome this limitation, we introduce a slight modification that ensures the discontinuity detector always considers every dimension even after the effective dimensionality of the manifold has been found. 3.3. Domain classification The grid generated by the discontinuity method can be used to classify the computational domain into its continuous subdomains. To obtain this decomposition the points in the grid must first be classified. The points are classified iteratively using a ray tracing algorithm. Starting with a randomly selected grid point the left and right neighboring points in each axial direction are found. If these points exist and are not separated by a discontinuity point they must belong to the same class as the original point. Any such connected neighbors are given the same classification as the original point and their connected neighbors identified. Once all points connected to the initial point have been found another point not yet classified is chosen and used as the basis for another classification cycle. The cost of domain decomposition method is Oðd nÞ where n is the number of points generated by the detector. Each point is visited only once, but each time a point is visited every dimension must be checked for neighbouring points. This approach is inexpensive but lacks robustness and a better classification algorithm would further enhance the utility of the minimal multi-element method. For the ray tracing algorithm to work the discontinuity detector must correctly identify a discontinuity between all points in the grid it generates. If it fails to do so the ray tracing algorithm will incorrectly combine two disjoint regions into one element. Constructing a classification algorithm that does not suffer from this deficiency is the subject of future work. Once all the points in the grid have been classified, any number of existing classification techniques can be used to construct a classification function to make predictions for new, yet unknown, data. In this paper we use a simple heuristic to

796

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

classify arbitrarily positioned points. The class of a new point is simply the class of the nearest grid point. This will incorrectly classify points in the d neighborhood of the discontinuity but this is within the resolution level specified by the user. Given a resolution d this classification procedure has been shown to provide an accurate classification function [20]. The accuracy of the classification and distance functions is proportional to the resolution used to resolve the discontinuity. The number of grid points is proportional to the dimension and hyper-dimensional surface area of the functional edge and inversely proportional to the resolution used to resolve the discontinuity. When the dimensionality of the discontinuity resides in a small subset of dimensions of the input space the proposed method is very efficient. When the dimensionality of the manifold is small, and the function is linear in the variables in the dimensions not belonging to the manifold, the algorithm becomes ‘‘optimal,’’ that is the total number of points required grows linearly with the dimensionality [20]. 4. Minimal-element method: element-level approximation In our mE method, since the discontinuity determines the domain decomposition, the resulting elements can be arbitrarily shaped, and thus traditional approximation techniques on cardinal domains such as hyper-cube are no longer appropriate. Though in principle these irregular elements could be mapped to a cardinal domain, the transformations to achieve this can be extremely complicated and challenging to implement, especially in high dimensions. Here we propose to use stochastic collocation (SC) to construct the element-level approximation. In particular, since our discontinuity detection algorithm already produces solution simulation results on a set of adaptive sparse grids, we will utilize this solution ensemble and will not perform additional simulations. The difficulty is that any structure of the underlying sparse grids, or any other structured collocation grids, is lost because of the irregular shape of the elements. Therefore we need a method to construct accurate gPC approximations using arbitrary number of collocation nodes located on arbitrary locations. The goal can be described in the following abstract manner: Given a set of distinct points fxðjÞ gNj¼1 , arbitrarily distributed in d-dimensional Euclidean space, and ffj ¼ f ðxðjÞ ÞgNj¼1 are samples of a target function f : Rd ! R, find a polynomial qðxÞ in a such that q  f in a well-defined sense. Here we propose to use interpolation for the approximation. By doing so we enforce the condition qðxðjÞ Þ ¼ fj ; j ¼ 1; . . . ; N. For arbitrary nodes in high dimensions, polynomial interpolation is not a trivial task. To this end, we employ the ‘‘least orthogonal interpolation’’ method developed in [25], largely because of many of its desirable mathematical properties. We remark that one is certainly free to choose other approximation methods. We have compared our proposed method with least-squares approximations and ‘1 -optimized sparse expansions and found that our method is generally superior. 4.1. The least orthogonal interpolant The material in this section is based upon the work contained in [25]. We include only a brief summary to illustrate the important concepts and properties of the method. We recall that our setup is N nodes X ¼ fxðnÞ gNn¼1 and function evaluations distributed on some element E in Rd with general geometry. 4.1.1. The interpolation space On the parameter element E , let the joint probability density function pðxÞ be strictly positive on the interior of E , with finite polynomial moments of all orders. On the space of functions square-integrable against the weight function p, let fUn g1 be a(ny) family of orthogonal polynomials, normalized to be orthonormal under the p norm: hUn ; Un ip ¼ R 2n¼1 Un pðxÞdx ¼ 1. We alternate between indexing the polynomials by natural numbers, as done above, and by d-dimensional E multi-indices a; the multi-index ‘1 norm kak corresponds to the total degree of Ua . Let P be the space of all d-variate polynomials. Pk is the space of polynomials of degree less than or equal to k. For any f 2 P, define the degree-k projection operator P k as

f ¼

X

X

f a Ua () P k f ¼

a

fa Ua ;

kak6k

where fa ¼ hf ; Ua ip . The goal is, given the N nodes X and the orthgonal polynomials generated from the density pðxÞ, to construct a space of polynomials PX;p with which we can uniquely perform interpolation from data on X. To accomplish this, we use the nodal locations to define a collection of generalized functions:

hn ðÞ ¼

X

Ua ðxðnÞ ÞUa ðÞ;

n ¼ 1; . . . ; N:

a

The algebraic functional kn ðÞ ¼ h; hn ip is point evaluation at xðnÞ . We introduce an operator that maps generalized functions h to polynomials:

h#;p ¼ P m h;

m ¼ minfk 2 N0 : P k h – 0g:

h#;p is a polynomial of (total) degree m. The explicit mention of the probability density p when defining the above operator indicates that indeed the choice of p matters since it determines the projection operators P k . We extend the operator ðÞ#;p to vector spaces H by H#;p ¼ spanfh#;p : h 2 Hg. If H ¼ spanfhn gNn¼1 one can show that the space of polynomials

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

797

PX;p ¼ H#;p is a dimension-N space of d-variate polynomials that is isomorphic to interpolation data at fxðnÞ g. Therefore, a polynomial g 2 PX;p satisfying gðxðnÞ Þ ¼ fn exists and is unique, and this polynomial is called the least orthogonal interpolant [25]. 4.1.2. Construction of the least orthogonal interpolant The construction of both the interpolation space and the necessary interpolation coefficients can be accomplished with standard numerical linear algebra operations on the interpolatory Vandermonde matrix. Let V be a matrix with rows indexed by natural numbers and columns indexed by multi-indices, with the columns ordered by total degree.4 V has entries V n;a ¼ Ua ðxðnÞ Þ. The algorithm originally presented in [9], and extended in [25] performs LU-type operations on V to obtain the decomposition

V ¼ LUH; where L and U are size-N square, invertible matrices, and H has the same size as V. As the variable names suggest, L and U are lower- and upper-triangular, respectively. The matrix H has a block-diagonal structure with orthogonal rows and contains transformation coefficients connecting the Un to a polynomial basis for PX;p . The interpolation coefficients in the basis on PX;p are given by c ¼ U1 L1 f, where f is the vector of data evaluations. Information from H can be used to easily transform c into coefficients for Ua . Because we cannot a priori predict the polynomial degree required for interpolation, the number of columns for V that are necessary is not known before the algorithm is started. However, one needs only to construct small submatrices of the dense matrix V at any particular stage; furthermore, these submatrices need only be generated once. Thus, even for high dimensional, high degree situations, the algorithm is feasible. 4.2. An adaptive interpolation scheme In the general setup of this paper, we will frequently have far more interpolation nodes (obtained from the discontinuity detector) than we need, and choosing a subset of them for interpolation will produce satisfactory results, assuming some smoothness of the function. The motivation for using a subset of the points is as follows: the discontinuity detector does not choose points that are useful for interpolation – rather, points are chosen in an effort to classify regions. Thus, attempting full interpolation on a ‘bad’ set of interpolation nodes will produce an ill-behaved result. We choose a subset of the given nodes that produces a well-behaved and informative interpolant. The methodology for doing is explained below. One attractive property of the least orthogonal interpolant scheme is that the polynomial spaces that are constructed are e X be another collection of nodes that contains the nodes from X. Then P PX;p . Therefore, nested. In other words, let X e X ;p once the work has been done to construct the interpolation space for X, we can ‘reuse’ the results to generate the interpoe e lation space for X . This can result in a substantial amount of savings, especially if X n X has low cardinality. We can use this property to form a computationally feasible method for interpolation that depends nonlinearly on the data. Let X be some proper subset of fxðnÞ g. To begin, we remark that basis coefficients for PX;p given by c ¼ U1 L1 f correspond to basis elements that are themselves polynomials orthonormal under h; ip . Because of this, the ‘2 norm of the coef: ficients kck (with c1 set to 0) is the standard deviation of the interpolant. Then for any i 2 I¼fxðnÞ g n X inclusion of this point into the interpolation algorithm results in the step:

e e 1 d; e c¼U e¼e d L 1 f; where

e L¼



L

0

0 1



INN ~lT

0



; e L Nþ1;Nþ1

e differs from U only in its last column. Then the first N entries with ~lT the appropriate submatrix of e L. Similarly, the matrix U e are identical to the entries for d at step N and need not be recomputed. Only the back-substitution application of U e 1 of d must be recomputed. In this way, we can compute candidate e c vectors for each i 2 I, and we choose the next node for interpolation as the i that minimizes the ‘2 norm of c (where c1 is set to zero).5 Though more computationally intensive than a straightforward application of the least interpolant method, this method is still very inexpensive. Even if the choice of new interpolation nodes begins to dominate computations, the part of the algorithm that chooses the next node is ‘‘embarrassingly’’ parallel. 4 An intra-degree ordering for multi-indices is not relevant for any of the algorithms, nor for our presentation. One may choose whatever ordering is convenient, so long as it adheres to a non-decreasing total degree ordering. 5 Naturally one may instead choose to minimize the ‘1 norm of the coefficients instead of the ‘2 norm.

798

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

Fig. 3. Two-dimensional surfaces of the four benchmark problems.

4.3. Minimal-element method: summary We are now in a position to finalize the minimal-element method we propose: first we use the discontinuity classifier described in Section 3 on the unknown functions f to decompose the domain into N e elements, with the kth element denoted E k . This requires only collocation evaluations of the unknown. Let Q k ; k ¼ 1; 2; . . . ; N e , denote the number of collocation evaluations of f that are enclosed by each element. For each element E k , we reuse the Q k collocation evaluations to reconstruct the polynomial surrogate fE k using the adaptive least orthogonal interpolant from Section 4.2. This defines the interpolation operator IE k in (2). Thus construction of fE k is a postprocessing step. The full piecewise-defined surrogate is then defined by (2). 5. Numerical examples In this section we present multiple examples that illustrate the effectiveness of the proposed minimal-element approximation collocation method. In the examples presented we consider a discontinuous function f : C ! R with d input variables x ¼ ðx1 ; . . . ; xd Þ that are defined on the input space C ¼ ½1; 1d with uniform weight p  1. The construction of orthogonal polynomials on a general geometry E k is a difficult problem itself and is not the topic of this paper. To compute the least orthogonal interpolant on E k , we construct a hypercube in Rd that circumscribes the Q k collocation nodes and use the orthogonal polynomials associated with that hypercube with uniform weight (e.g. tensor-product Legendre polynomials).

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

799

Fig. 4. Collocation nodes (black) generated by the discontinuity detector. Red points represent intervals identified as containing a discontinuity. For the cube, sphere, and triagle, there are Q 1 points used on the interior, and Q 2 points on the exterior. For the plane, there are Q 1 points on the left-hand side, and Q 2 points on the right-hand side. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

To determine the accuracy between f and fEk , we use Monte Carlo estimates. First we generate 10; 000 Monte Carlo samples on ½1; 1d . The discontinuity detector classifies each of these parameter samples into the appropriate region. Some samMonte Carlo points ples are discarded when they are within a small distance from the discontinuity.6 We thus have Nerror Ek Nerror k

E for the element E k , and the ‘1 , ‘2 , and ‘1 errors for each region are then approximated by fxE k ;n gn¼1

kf  fEk k1 ’ kf  fEk k22 ’

1

Nerror k

     E X   f xE k ;n  fEk xE k ;n ;

Nerror Ek n¼1 1

Nerror k

    2 E X   f xEk ;n  fE k xE k ;n  ;

Nerror Ek n¼1

       kf  fEk k1 ’ maxf xE k ;n  fEk xE k ;n : n

6 Discarding points in such cases is motivated by the facts that the distance of these points to the discontinuity manifold lies below the d resolution of the discontinuity detector, which calls into question the fidelity of the classification

800

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

Remark 1. The multi-element method presented here has been formulated to approximate functions with jump discontinuities in a function. If one wishes to approximate functions with discontinuities in derivatives the polynomial annihilation method presented in Section 3.1 must be reformulated as discussed in [4]. The adaptive algorithm which facilitates the application of discontinuity detection to high dimensional applications and the interpolation procedure,

Fig. 5. Classification of Monte-Carlo samples by the discontinuity detector for the four benchmark problems. Red points cannot be classified with certainty as they fall within the d resolution employed by the detector. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6. Errors for the minimal element method on test function f cube .Hereafter,

DOF

indicates the number of samples used for interpolation.

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

801

however, do not need to be changed. The local stencils generated by the adaptive algorithm, and employed by the polynomial annihilation method (3), can still be used, only the coefficients cj ðxÞ and normalization factor qm ðxÞ in must be recalculated.

Fig. 7. Errors for the minimal element method on test function f sphere .

Fig. 8. Errors for the minimal element method on test function f triangle .

Fig. 9. Errors for the minimal element method on test function f2plane .

802

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

Fig. 10. Contour plots for the minimal element approximation. Each pair of contour plots shows (top) the approximation with the indicated number of points and (bottom) the exact contour lines with interpolation points chosen for the adaptive least orthogonal method indicated by dots. N 1 and N 2 indicate the number of points used for interpolation in regions 1 and 2, respectively.

5.1. Benchmark problems To illustrate the numerical properties of the proposed multi-element method we consider standard benchmark problems from the Genz function set [16], with discontinuities imposed on certain manifolds. Fig. 3 shows the surfaces of the twodimensional versions of these benchmark functions.

( f cube ðxÞ ¼

1 þ f 1 ðxÞ; kx1 k < 0:45 and kx2 k < 0:45; f1 ðxÞ;

( f sphere ðxÞ ¼

1 þ f 1 ðxÞ;

Pd

f 1 ðxÞ;

otherwise:

( f

triangle

ðxÞ ¼

f 1 ðxÞ;

ð7aÞ

otherwise: 2 i¼1 xi

< 0:52 ;

x1 þ x2 > 14

and x1  x2 > 14

ð7bÞ

and x1 < 0:5;

f 1 ðxÞ þ 3; otherwise:

Fig. 11. Errors for the minimal element method on test function f6plane .

ð7cÞ

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

( fdplane ðxÞ ¼

f 1 ðxÞ;

3x1 þ 2x2  0:01 > 0;

fd2 ðxÞ þ 12 cosðpðx1 þ x2 þ 0:3ÞÞ þ 1; otherwise:

803

ð7dÞ

where 2 X f ðxÞ ¼ exp  x2i 1

!  x31  x32 ;

i¼1

fd2 ðxÞ ¼ 1 þ f 1 ðxÞ þ

d 1 X x2 : 4d i¼2 i

Fig. 4 displays the collocation nodes generated by the discontinuity detector for each of the four two-dimensional benchmark problems. Red points represent intervals identified as containing a discontinuity. The classification of the Monte-Carlo samples for the four benchmark problems are visualized in Fig. 5. Note that the classification does not require functional evaluations. These two-dimensional functions are decomposed into two elements. Each sample is either classified as belonging to element E 1 (blue points) or element E 2 (green points). Red points represent samples that cannot be classified with certainty as they fall within the d resolution employed by the detector. The discontinuity detector has two main tunable parameters, the minimum resolution d and the maximum order of the polynomial annihilation scheme. In all the examples d ¼ 2  106 and the maximum order in (6) is m ¼ 6.

Fig. 12. Discontinuity detection applied to a two dimensional function f2multi with multiple discontinuities. The true function is shown in (a), the points generated by the discontinuity detection algorithm are shown in (b), and 10,000 randomly classified points are shown in (c). Four smooth regions are identified. Green points represent the subset of random points that cannot be classified with certainty. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

804

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

5.2. Two domains, two dimensions In this section we illustrate the effectiveness of the minimal element method for response surface reconstruction. Figs. 6– 9 show the errors for the four test functions in (7) when d ¼ 2. Using the discontinuity detector, the input domain ½1; 12 of each function is decomposed into two elements. For each function a high level of accuracy is achieved with only a very small subset of the available collocation nodes. When applied to each d ¼ 2 example from (7), the accuracy of the interpolant in both regions increases as the number of the collocation nodes used increases. This behavior observed is similar for all of the two-dimensional geometries. Fig. 10 shows how the minimal element approximation, when applied to (7d), changes as the number of collocation nodes increases. It also shows how the adaptive scheme chooses points for interpolation.

5.3. Comparison of interpolation methods Since the polynomial surrogate reconstruction procedure is a post-processing operation, one can envision using any of a number of interpolation or regression methods to construct such a surrogate. For the two-domain, two-dimensional triangular example (7c) we have run tests that compute (a) the ‘2 minimum-norm/least-squares solution, (b) the interpolant using ‘1 norm minimization (‘basis pursuit’) on the coefficient vector, (c) relaxed interpolation with ‘1 norm minimization (‘basis pursuit denoise’), and (d) a straightforward (non-adapted) implementation of the least orthogonal interpolant. Our results indicate that all methods other than the least-squares approach are less robust. However, when compared to the least-squares approach, the adaptive least orthogonal method has the advantage that the solution space is progressively built up degree-by-degree. This means that one only requires submatrices of the interpolatory Vandermonde matrix, rather than the full matrix. Additionally, Information from previous degrees is used to construct the interpolant for higher degrees [25]. In contrast, the least-squares method (a) requires information from the full interpolatory Vandermonde matrix, and (b) one must recompute the entire solution if the solution basis is changed.

Fig. 13. Errors for the minimal element method on test function f2multi .

805

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

5.4. Approximation in six dimensions Consider now fdplane with d ¼ 6. Again the discontinuity detector decomposes the input space into two elements. The adaptive least orthogonal method is then used to construct an interpolant on each element using the set of collocation nodes generated by the discontinuity detector. In this case running some of the alternatives enumerated in the previous section become prohibitively expensive: the space of fifth-degree sexavariate polynomials has dimension 426 and simply storing, or evaluating matrix-vector products with the interpolatory Vandermonde matrix on all Q k points becomes problematic.

Fig. 14. Errors for the minimal element method on test function f6multi .

Table 1 ‘2 error in the minimal element and sparse grid approximation on the two elements E 1 and E 2 for examples (7a)–(7d) and (8). N (left) is the total number of function evaluations used by the discontinuity detector, and the Minimal Element interpolation method uses a total of DOF evaluations for interpolation over all elements. N (right) is the number of points used by the adaptive sparse grid method, with the specified error tolerance. The ‘‘tol’’ factor is the tolerance level used for the adaptive sparse grid approximation. Function

N

Minimal element 2

N

1

2

2

DOF

‘ ðE Þ

‘ ðE Þ

Adaptive sparse grid tol

‘2 ðE 1 Þ

‘2 ðE 2 Þ

f sphere

655

200

1:52  109

3:59  105

6285

1:00  104

1:24  102

4:32  102

f cube

247

149

5:95  107

3:42  106

1297

1:00  104

4:35  104

3:25  103

f triangle

930

200

1:15  108

1:02  104

11,239

1:00  104

5:66  102

1:12  101

f2plane

829

200

1:69  105

1:42  104

7153

1:00  104

1:98  102

4:92  102

f6plane f6plane

7804

200

9:19  10

2

2

3

2

7804

300

1:75  102

1:96  10

45,601

1:00  10

2:22  10

5:39  102

4:77  103

225,443

1:00  104

2:17  103

2:43  102

806

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

In contrast, the 100-DOF adaptive least orthogonal interpolant is also a fifth-degree polynomial but its construction requires storage of only a few 100  100 matrices, with one-time construction and sequential storage of small submatrices of V. Fig. 11 shows results for the plane discontinuity using the adaptive least orthogonal method. As with the two-dimensional examples the error initially decreases until an optimal number of collocation nodes has been employed. After this point, the introduction of additional nodes slightly degrades the accuracy of the interpolant. In this manuscript the choice of how many points to use in the construction of the interpolant is arbitrary. Adding a modification to the algorithm that stops adding points when the accuracy of the interpolant decreases is the topic of future work.

5.5. Approximating functions with multiple discontinuities The minimal element collocation method proposed in this manuscript is not restricted to functions that possess only one discontinuity. It can also be used for functions with multiple discontinuities. Consider the function

8 > f 1 ðxÞ  2; > > > < 2f 2 ðxÞ; d fdmulti ðxÞ ¼ > 2f 1 ðxÞ þ 4; > > > : 1 f ðxÞ;

3x1 þ 2x2 P 0 and

 x1 þ 0:3x2 < 0;

3x1 þ 2x2 P 0 and

 x1 þ 0:3x2 P 0;

ðx1 þ 1Þ2 þ ðx2 þ 1Þ2 < 0:952

ð8Þ

and d ¼ 2;

otherwise:

The function surface is shown in Fig. 12(a). When applied to this function, the minimal element method splits the input domain ½1; 12 into four elements. The collocation nodes generated by the discontinuity detector are shown in Fig. 12(b) and the classification of 10; 000 random Monte Carlo samples is shown in Fig. 12(c). Fig. 13 displays the error in each of the four elements obtained when the adaptive least orthogonal method is used to construct the interpolants. In all regions a very high level of accuracy is achieved using only a small number of points. Similar results are obtained when the minimal element collocation method is applied to the six-dimensional problem f6multi ðxÞ (refer to Fig. 14), although the final accuracy is lower than in the two-dimensional case.

5.6. Comparison with adaptive sparse grid collocation In the previous sections we have shown the utility of the proposed minimal element collocation method for approximating discontinuous functions. Here we compare the proposed method to the popular adaptive sparse grid collocation method [22]. We note that there exist a few recent works on discontinuity detection in high dimensions, e.g. [11,26]. These are relatively new methods. Here we have chosen the adaptive sparse grid method for comparison because of its general applicability and popularity. Many of the discontinuity detection methods either have difficulty representing complicated manifolds, or cannot resolve multiple regions, and so direct comparison with our method is difficult. While the sparse grid method was not designed to handle discontinuities in particular, it performs relatively well in their presence. For comparison against the sparse grid method we choose the standard error criterion based upon the magnitude of the hierarchical surplus and the volume of the associated basis function. Table 1 shows the ‘2 error in the minimal element and sparse grid approximations on the two elements E 1 and E 2 of the four two-dimensional benchmark problems (7a)–(7d) and the six-dimensional example f6plane . Table 2 compares the error in the minimal element and sparse grid approximations when applied to the functions f2multi and f6multi with multiple discontinuities. In almost all cases the minimal element method achieves superior levels of accuracy compared to the sparse grid method, with significantly fewer points. Indeed in many cases there is at least an order of magnitude difference. The difference in the two methods efficiency is smallest for the function (7a). This can be directly related to the nature of the underlying discontinuity. In this function the discontinuity lies parallel to the axial directions, for which the sparse grid method is ideally suited. Despite this fact, the multi-element method still performs considerably better.

Table 2 ‘2 error in the four elements E 1 ; . . . ; E 4 for the two dimensional function f2multi with multiple discontinuities. Function

‘2 ðE 1 Þ

‘2 ðE 2 Þ

‘2 ðE 3 Þ

‘2 ðE 4 Þ 3:72  105 –

Minimal element 2D: 1522 points, 200 DOF. 6D: 9856 points, 300 DOF f2multi

3:69  107

7:12  104

9:39  104

f6multi

7:42  104

1:83  102

2:70  102

Adaptive sparse grid 2D: 19474 points, 6D: 320476 points f2multi

3:56  102

2:17  102

7:21  102

f6multi

7:03  102

2:13  102

3:43  102

5:52  102 –

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

807

6. Conclusion In this paper we present a multi-element interpolation method that can be used to approximate discontinuous multivariate functions. The method employs discontinuity detection to decompose the input space into disjoint regions of high-regularity. If any discontinuities are present, they are used to naturally decompose the domain into a minimal number of elements. Once the domain has been decomposed, we then construct a polynomial interpolant on each element. Because we form elements using locations of function discontinuities as boundaries, the shape of the regions constructed are often irregular. Consequently traditional approximation methods based upon structured sets of nodes are impractical. To overcome this difficulty we construct an approximation scheme that can construct a gPC approximation in each smooth region from an arbitrary distribution of interpolation nodes. This scheme removes the traditional restriction to hyper-cube elements of traditional multi-element methods. Moreover it allows the use of the nodes generated by the discontinuity detection algorithm to be reused to construct the local interpolant. The interpolant on each element is built using an adaptive variation of the least orthogonal interpolant. The size of the polynomial basis produced is equal to the number of points used to construct the interpolant. To increase the efficiency of the method we build the polynomial approximation on a subset of the points created by the discontinuity detector. This choice of subset is driven by a greedy variance-suppressing algorithm that is efficient due to the nested basis property of the least orthogonal interpolant. We apply the proposed minimal element approximation method to a number of multivariate discontinuous functions and show that regardless of the shape or number of discontinuities the method achieves high levels of accuracy with only a small number of function evaluations. When compared to the more established adaptive sparse grid collocation scheme, the proposed method achieves a higher order of accuracy with significantly fewer function evaluations. In some cases there is an order of magnitude difference in the number of points used by the two methods.

References [1] N. Agarwal, N.R. Aluru, A domain adaptive stochastic collocation approach for analysis of mems under uncertainties, J. Comput. Phys. 228 (20) (2009) 7662–7688. [2] R. Archibald, A. Gelb, R. Saxena, D. Xiu, Discontinuity detection in multivariate space for stochastic simulations, J. Comput. Phys. 228 (7) (2009) 2676– 2689. [3] R. Archibald, J. Gelb, A. Yoon, Polynomial fitting for edge detection in irregularly sampled signals and images, SIAM J. Numer. Anal. 43 (1) (2005) 259– 279. [4] R. Archibald, J. Gelb, A. Yoon, Determining the locations of discontinuities in the derivatives of functions, Appl. Numer. Math. 58 (5) (2008) 577–592. [5] I.M. 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. [6] I.M. Babuska, R. Tempone, G.E. Zouraris, Galerkin finite element approximations of stochastic elliptic partial differential equations, SIAM J. Numer. Anal. 42 (2) (2004) 800–825. [7] R. Bauer, Band Pass Filters for Determining Shock Locations, PhD thesis, Applied Mathematics, Brown University, 1995. [8] M. Bieri, C. Schwab, Sparse high order FEM for elliptic SPDEs, Comput. Methods Appl. Mech. Eng. 198 (2009) 1149–1170. [9] C. De Boor, A. Ron, Computational aspects of polynomial interpolation in several variables. Mathematics of Computation, 58(198):705–727, April 1992. ArticleType: primary_article / Full publication date: Apr., 1992 / Copyright –1992 American Mathematical Society. [10] J. Canny, A computational approach to edge detection, IEEE Trans. Pattern Anal. Machine Intell. 8 (1986) 679–698. [11] T. Chantrasmi, A. Doostan, G. Iaccarino, Padé-legendre approximants for uncertainty analysis with discontinuous response surfaces, J. Comput. Phys. 228 (19) (2009) 7159G–7180. [12] J. Foo, G.E. Karniadakis, Multi-element probabilistic collocation method in high dimensions, J. Comput. Phys. 229 (5) (2010) 1536–1557. [13] J. Foo, X. Wan, G.E. Karniadakis, The multi-element probabilistic collocation method (ME-PCM): error analysis and applications, J. Comput. Phys. 227 (22) (2008) 9572–9595. [14] T. Gardner, C. Cantor, J. Collins, Construction of a genetic toggle switch in escherichia coli, Nature 403 (2000) 339–342. [15] A. Gelb, E. Tadmor, Adaptive edge detectors for piecewise smooth data based on the minmod limiter, J. Sci. Comput. 28 (2-3) (2006) 279–306. [16] A. Genz, Testing multidimensional integration routines, in: Proceedings of International Conference on Tools, Methods and Languages for Scientific and Engineering Computation, Elsevier North-Holland, Inc., New York, NY, USA, 1984, p. 81G–94. [17] R.G. Ghanem, P.D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer-Verlag, New York, Inc., New York, NY, USA, 1991. [18] M. Griebel, Adaptive sparse grid multilevel methods for elliptic PDEs based on finite differences, Computing 61 (2) (1998) 151–179. [19] J.D. Jakeman, Numerical Methods for the Quantification of Uncertainty in Discontinuous Functions of High Dimension. Ph.D thesis, The Australian National University, November 2010. [20] J. D Jakeman, R. Archibald, D. Xiu, Characterization of discontinuities in high-dimensional stochastic problems on adaptive sparse grids, J. Comput. Phys. 230 (10) (2011) 3977–3997. [21] J.D. Jakeman, S.G. Roberts, Local and dimension adaptive stochastic collocation for uncertainty quantification, in: Jochen Garcke, Michael Griebel (Eds.), Sparse Grids and Applications, Lecture Notes in Computational Science and Engineering, vol. 88, Springer, Berlin Heidelberg, 2013, pp. 181–203. [22] 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. [23] X. Ma, N. Zabaras, An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations, J. Comput. Phys. 229 (10) (2010) 3884–3915. [24] L. Mathelin, M. Hussaini, T. Zang, Stochastic approaches to uncertainty quantification in CFD simulations, Numer. Algorithms 38 (1-3) (2005) 209–236. [25] A. Narayan, D. Xiu, Stochastic collocation methods on unstructured grids in high dimensions via interpolation, SIAM J. Sci. Comput. 34 (3) (2012) A1729–A1752. [26] K. Sargsyan, C. Safta, B. Debusschere, H. Najm, Uncertainty quantification given discontinuous model response and a limited number of model runs, SIAM J. Sci. Comput. 34 (1) (2012) B44–B64. [27] I. Sobel, An isotropic 3 image gradient operator, in: H. Freeman (Ed.), Machine Vision for Three-Dimensional Scenes, Academic Press, Boston, 1990. [28] M. Tatang, W. Pan, R. Prinn, G. McRae, An efficient method for parametric uncertainty analysis of numerical geophysical model, J. Geophys. Res. 102 (D18) (1997) 21925–21932.

808

J.D. Jakeman et al. / Journal of Computational Physics 242 (2013) 790–808

[29] X. Wan, G.E. Karniadakis, An adaptive multi-element generalized Polynomial Chaos method for stochastic differential equations, J. Comput. Phys. 209 (2) (2005) 617–642. [30] X. Wan, G.E. Karniadakis, Multi-element generalized Polynomial Chaos for arbitrary probability measures, SIAM J. Sci. Comput. 28 (3) (2006) 901–928. [31] D. Xiu, Numerical Methods for Stochastic Computations: A Spectral Method Approach, Princeton University Press, 2010. [32] D. Xiu, J.S. Hesthaven, High-order collocation methods for differential equations with random inputs, SIAM J. Sci. Comput. 27 (3) (2005) 1118–1139. [33] D. Xiu, G.E. Karniadakis, The Wiener–Askey Polynomial Chaos for stochastic differential equations, SIAM J. Sci. Comput. 24 (2) (2002) 619–644.