An Adaptive WENO Collocation Method for Differential Equations with Random Coefficients

The stochastic collocation method for solving differential equations with random inputs has gained lots of popularity in many applications, since such a scheme exhibits exponential convergence with smooth solutions in the random space. However, in some circumstance the solutions do not fulfill the smoothness requirement; thus a direct application of the method will cause poor performance and slow convergence rate due to the well known Gibbs phenomenon. To address the issue, we propose an adaptive high-order multi-element stochastic collocation scheme by incorporating a WENO (Weighted Essentially non-oscillatory) interpolation procedure and an adaptive mesh refinement (AMR) strategy. The proposed multi-element stochastic collocation scheme requires only repetitive runs of an existing deterministic solver at each interpolation point, similar to the Monte Carlo method. Furthermore, the scheme takes advantage of robustness and the high-order nature of the WENO interpolation procedure, and efficacy and efficiency of the AMR strategy. When the proposed scheme is applied to stochastic problems with non-smooth solutions, the Gibbs phenomenon is mitigated by the WENO methodology in the random space, and the errors around discontinuities in the stochastic space are significantly reduced by the AMR strategy. The numerical experiments for some benchmark stochastic problems, such as the Kraichnan-Orszag problem and Burgers’ equation with random initial conditions, demonstrate the reliability, efficiency and efficacy of the proposed scheme.


Introduction
Problems subject to uncertainty arise in many engineering [1,2], environmental and biological applications.Such uncertainty is mainly due to a lack of knowledge about the true value of model parameters or the random nature of the quantity of interest being studied.Uncertainty quantification (UQ) for practical problems has been drawing growing interest in recent years, in particular in developing numerical methods for stochastic computations.Due to the computational cost of high-fidelity models for complex systems, many researchers are investigating more efficient stochastic algorithms than the classic Monte Carlo (MC) method.One of the popular stochastic methodologies is generalized polynomial chaos (gPC) [3], which is an extension of the standard polynomial chaos method [4].In addition, advanced gPC algorithms for high-dimensional stochastic problems have been developed [5][6][7][8][9][10][11][12][13][14][15].
The performance and convergence of gPC depends on the smoothness of the function in the random parameter space.gPC will converge quickly for smooth functioning in random parameter space.However, the performance and convergence will deteriorate when the function has low regularity, particularly with discontinuities in the random parameter space.In such cases, Gibbs-type phenomenon will be observed in gPC approximation, which causes the slow convergence.To alleviate such a challenging issue, multi-element gPC methods [16,17] and multi-element stochastic collocation methods [18,19] were developed to decompose the random space into sub-domains, then employ a gPC expansion in each element.The challenge of the standard multi-element gPC or multi-element stochastic collocation method is to detect the discontinuities in high-dimensional random space and the computational cost.It requires the users to resolve the stochastic problem in each of the random element.A few recent works on the development of efficient stochastic algorithms for handling discontinuities also exist [20][21][22][23][24][25][26].
In this paper, we propose an adaptive high-order multi-element stochastic collocation method in conjunction with Weighted Essentially non-oscillatory (WENO) reconstruction methodology and adaptive mesh refinement (AMR) strategy to address the challenging issue regarding standard multi-element approaches for solving stochastic differential equations with discontinuities in the random space.The proposed method consists of two important components.(1) Adaptive stochastic domain decomposition framework: it is known that solution in the random space may exhibit local complex structures such as discontinuities; hence, we employ a stochastic domain decomposition framework based on the idea that in the smooth region a multi-element stochastic collocation scheme is used in order to take advantage of its low computational cost and fast convergence, while an AMR WENO scheme on a local tensor-product grid is used in the region with complex discontinuous solution structures; (2) High-order WENO interpolation scheme in random space: it is a grand challenge to detect the discontinuities in high-dimensional random space.However, by employing the WENO scheme, there is no need to know the location of the discontinuities exactly.High-order WENO interpolation scheme in random space can provide essentially non-oscillatory solutions across the discontinuities, which is particular suitable for stochastic problems with complex discontinuity structure; (3) AMR strategy in random space [27]: Even though the AMR WENO scheme has been developed in different settings in the literature, such as an AMR WENO scheme for solving hyperbolic conservation laws [27], it has not been investigated in the uncertainty quantification context.By employing AMR strategy, we can identify the regions that contain discontinuities in random space and only apply the high-order WENO interpolation scheme in such regions.In this paper, we will show the efficiency and efficacy of the proposed scheme for solving stochastic problems with complex solution structure including discontinuities.
This paper is organized as follows: in the next section, the mathematical formulation of the proposed stochastic domain decomposition framework, high-order WENO interpolation, and AMR strategy in random space is introduced.In Section 3, the numerical examples are given.Finally, concluding remarks are provided in Section 4.

Formulation
In this section, we formulate the adaptive AMR WENO-based multi-element stochastic collocation scheme for solving differential equations with random inputs.We start with a review of the multi-element stochastic collocation method, the high-order WENO interpolation methodology for non-smooth problems, which has been widely used in many applications.In the second part, an ARM strategy in conjunction with the WENO interpolation is formulated in order to take advantage of both numerical ingredients.

Multi-Element Stochastic Collocation Method (ME-SCM)
The ME-SCM introduced in [19] decomposes the random space Γ into non-overlapping elements and then employs the stochastic collocation method on each element.This will provide the approximation of local moment statistics in each element.We can assemble such local moment statistics in each element to obtain the global statistics.To improve the efficiency of the proposed ME-SCM, both the high-order WENO interpolation scheme and AMR strategy will be employed, which will be discussed in the following subsections.
In particular, we discretize the parametric space Γ into a nonoverlapping mesh of open hypercubes.We denote {B i } N e i=1 to be a finite collection of open subsets of Γ such that N e i=1 Bi = Γ and B i B j = ∅ whenever i = j.For illustrative purposes, we let B i be hypercubes.These sets consists of the elements of a mesh on the random space Γ and N e denotes the number of elements.In the prescribed mesh, a set of collocation points {y i,q } s q=1 is prescribed in each element B i .Here, s refers to the number of collocation points.The collocation points are selected so that they coincide with the points of an cubature rule on B i with integration weights {w i,q } s q=1 .In particular, we employ full tensor products of Gauss quadrature points as the collocation points.
At each of the collocation points {y i,q } s q=1 , we obtain the numerical solution u k (x, y i,q ) of the corresponding deterministic problem, L(x, y i,q ; u k (x, The fully discrete approximation I B i u k (x, y) based on the tensor product Lagrangian interpolant can be written as: where p determines the degree of the interpolant in each dimension.l i j (y) is the Lagrange polynomial with respect to the quadrature point y i,q .
The global approximant is defined as, where Υ {y∈B i } is the characteristic function of set B i .The corresponding probability density function in each element is defined as, where ρ(y) = ∏ N j=1 ρ j (y j ).The local mean of a function u in an element i is given by, The approximation of local mean of ũ using the cubature rule over each element can be computed as, where Ẽi is the expected value approximation operator using numerical quadrature rule.
The approximate global mean can be obtained from the local mean as: Other statistics can be obtained by a similar procedure.

High-Order WENO Interpolation in Random Space
The widely recognized WENO interpolation method is known to be able to effectively resolve non-smooth data.By carefully assigning weights to each stencil candidate, the WENO interpolation attains high order accuracy in the smooth region while avoiding oscillations in regions with singularities (see [28] for a detailed review on the recent development and applications).

One-Dimensional Case
We start with the simple one-dimensional case.The extension to high-dimensional problems can be established by a tensor product construction and will be discussed later.Consider the following problem: given a set of point values u i on a uniform mesh x i in a bound domain, i = 1, . . ., N with a mesh size ∆x = x i+1 − x i , we need to interpolate a local polynomial approximating the original function u(x) with high order accuracy and free of oscillations.Such a local polynomial u h i (x) can be reconstructed as follows.
Assume that we would like to obtain a fifth order approximation from five point values u i−2 , u i−1 , u i , u i+1 , u i+2 .Rather than directly reconstructing a polynomial on the interval [x i−1/2 , x x+1/2 ], we first interpolate the values at five local quadrature points x i,q , q = 1, . . ., 5, see Figure 1, then obtain the polynomial from the values at these points.Such an idea has been used to develop the WENO limiter for discontinuous Galerkin (DG) schemes [29].First, note that for each point x i,q , we have three interpolation values 2,q u i+2 from three small stencils as third order approximations, and v i,q = c 0,q u i−2 + c 1,q u i−1 + c 2,q u i + c 3,q u i+1 + c 4,q u i+1 + c 5,q u i+2 from the large stencil as a fifth order approximation of the underlying function.We can also seek a set of linear weights d (r) q , r = 0, 1, 2, such that Note that a direct application of linear weights d (r) q will lead to oscillations due to the famous Gibbs phenomenon.Instead, they are replaced by the nonlinear weights, which is known as the main idea of the WENO methodology.Denote three local polynomials of degree two on each small stencil by p (r) (x).Then, the smoothness indicators measuring the relative smoothness of each stencil are defined as: In particular, it gives Then, the nonlinear weights are defined accordingly as: where ε is chosen as a small number 10 −6 in order to avoid zero denominator.Lastly, the interpolation values at quadrature nodes are obtained by applying the nonlinear weights, i.e., ṽi,q = v (0) i,q w (0) i,q w (1) i,q w (2) q (9) Based on the WENO interpolation values at the quadrature points ṽi,q , q = 1, . . ., 5, we are able to reconstruct a polynomial ũ(x) of degree four.We remark that ṽi,q , q = 1, . . ., 5 share the same set of smoothness indicators, and the interpolation polynomial ũ(x) automatically attains the desired essentially non-oscillatory properties.If a cell is close to boundaries, then we need to apply a one-sided WENO interpolation procedure proposed in [30] to avoid oscillations.For brevity, we do not review the details but refer the readers to [30].We also remark that, if only the first several moments such as mean or variance are of interest, the WENO reconstruction can be directly applied to compute such quantities to save computational cost.In particular, we only need to perform the WENO interpolation for several moments rather than the whole solution.

Two-Dimensional Extension
The extension of the WENO interpolation methodology to two-dimensional cases is based on a tensor product construction.Specifically, in order to obtain a fifth order approximation ũij (x), we need the stencil plotted in Figure 2. Similar to a one-dimensional case, we need to interpolate the point values ṽij,qg , q, g = 1, . . ., 5 on the local cell [x i−1/2 , x i+1/2 ] × [y j−1/2 , y j+1/2 ], then reconstruct a polynomial of degree 4. The interpolation procedure is performed in a dimension-by-dimension fashion.To obtain the value ṽij,qg (the red circle in Figure 2), we first apply the one-dimensional WENO interpolation in x-direction to get v ir,q , r = j − 2, . . ., j + 2 (the blue circles in Figure 2).Then, we can further apply the interpolation in y-direction.Once the values ṽij,qg are obtained, the interpolating polynomial is constructed from ṽij,qg .

AMR Methodology in Random Space
The AMR methodology is known as an efficient numerical tool for solving problems with local complex solution structures, which has soared in popularity in many applications since the 1980s [31,32].In this work, we consider a cluster based method to generate fine grids where the local mesh refinement is needed.It is also known that such a method is very suited for WENO interpolation on a tensor product grid.The AMR interpolation methodology has the following steps: 1. Use a multi-resolution analysis to flag the regions of refinement interest (RRI): We compute the difference between the solution u i,j and the average of its four neighbors b = 1 4 (u i+1,j + u i−1,j + u i,j+1 , u i,j−1 ).If |u i,j − b| > ∆x, where is a user specified constant, then the point and its eight neighbors are marked as regions of refinement interest (RRI).2. Generate finer grids using the algorithm proposed in [33]: We generate some non-overlapping rectangular clusters that cover the RRI.In each cluster, the grid is refined with a ratio of three, i.e., ∆x l+1 = ∆x l /3, where ∆x l denotes the mesh size on the mesh level l (see Figure 3; two non-overlapping clusters are generated, and red circles denote point values on the fine grid).3. Repeat steps 1 and 2 until the maximal mesh level is attained.4. Generate the WENO interpolating polynomial by applying the algorithm reviewed in Section 2.2.
We start from the clusters with the finest mesh level.The needed boundary values in the WENO interpolation procedure can be either obtained from the neighboring cluster with same mesh level or interpolated from coarser clusters.Again, the WENO interpolation should be used in order to avoid the oscillations (see Figure 3; the green circles are interpolated from the blue circles on a coarser cluster).Note that the refinement ratio is set as three, since the solution on the coarse mesh level can be reused.Even though we only consider the two-dimensional case, the algorithm in [33] can be applied to problem in arbitrary dimension.Further, the rectangular clusters may contain some nodes are not RRI, however, the cutoff ratio is set as 0.7 meaning that at least 70% nodes included in the clusters are RRI.
Lastly, we remark that the use of tensor product grid seems less efficient than the sparse grid for high-dimensional problems.However, one should note that for many real applications in UQ, such as discontinuous response, the solution may have local dense structures in the random space.An adaptive sparse-grid based stochastic collocation scheme such as the one proposed in [26] still needs sufficient nodes clustered in the region to adequately represent the solution structures.Hence, in such a scenario, the computation cost of the tensor product grid is comparable to that of the sparse grid.Another rationale of the tensor product grid is that it is quite complicated to implement a high order WENO interpolation on an adaptive sparse grid, which is highly non-uniform.We leave the investigation of the WENO interpolation in conjunction with a sparse grid to the future work.

Approximation Investigation for Non-Smooth Functions
We first assess the approximation properties of the proposed AMR WENO stochastic collocation scheme.Below, we consider two functions with singularities in domain [0, 1] 2 : 2 + cos(2π(x + y)) otherwise Note that f 1 has discontinuous derivative on the circle x 2 + y 2 = 0.3 and f 2 is discontinuous with a line singularity parallel to the grid line.We start with the a tensor-product mesh 50 × 50 and local mesh refinement is performed based on the multi-resolution analysis.The maximum mesh level is set as 4 and 5 for f 1 and f 2 , respectively.In Table 1, we report the "error in variance", that is, we compute the "numerically computed variance" f 2 (x, y)dxdy − ( f (x, y)dxdy) 2 using the adaptive WENO collocation method and compare the result with the "exact variance", and the maximum error at 1000 randomly sampled points with different for f 1 .It can be observed that the errors decrease as is getting smaller, while the number of sampling points is getting larger accordingly.Furthermore, if the mesh level is increased, the error is reduced by the local mesh refinement strategy.However, we remark that for a fixed , there is no gain in accuracy by excessively increasing the mesh level, since the error eventually will be dominated by the coarsest mesh.For instance, the error by the three-level mesh is comparable to that by the four-level mesh (see Table 1).Since the singular curve is 'dense' in the one-dimensional space, the sparse grid algorithm will generate a 'local' tensor product grid around the curve in order to resolve such a singularity of the solution.In Figure 4, we report the interpolant functions with the AMR WENO reconstruction method and the corresponding AMR meshes for = 10, 1.It is observed that the local mesh refinement is incurred around the line singularity of the function.Furthermore, the local refined region with smaller is greater than that with larger .We would like to use function f 2 to demonstrate the advantages of using WENO reconstruction Equation ( 9) over the linear reconstruction Equation ( 8), as shown in Figure 5.Note that the maximum error for such a discontinuous function remains O(1), and hence we only report the convergence study for variance.In Table 2, we summarize the error convergence with = 1 for both reconstruction methodologies.Similar to the previous example, as we increase the mesh level, the error is reduced accordingly for both methods.Moreover, the errors computed by the WENO reconstruction are about one half of errors by the linear reconstruction, which justifies the usage of the robust WENO reconstruction methodology.We also remark that the methods with = 0.1, 1, 10 will generate the same AMR mesh.This is because f 2 itself is discontinuous, and hence the local mesh refinement is only carried out around discontinuities.

Two-Dimensional Kraichnan-Orszag (K-O) Problem
In this subsection, we consider the following transformed two-dimensional K-O three-mode problem: subject to the random initial condition Similar to [26], we let Y 1 and Y 2 have the uniform distribution Y 1 , Y 2 ∼ U(−1, 1).The stochastic simulation for this benchmark problem is challenging since the solution exhibits a bifurcation on the parameter y 1 (0) and y 2 (0) [16].Note that a direct application of gPC will result in failure in convergence after a short period of time.Here, we would like to use this example to further demonstrate the performance of the AMR WENO stochastic collocation scheme.The maximum mesh level is set to four and a classic forth order Runge-Kutta method with time step 0.01 is used to solve the differential equations.First, we plot the numerical solution of y 1 in the random space at T = 10 with = 2, 10 and the associated adaptive mesh in Figure 6.The number of points sampled in the simulation is 29,460 and 3684 for = 2, 10, respectively.It is observed that the local mesh refinement is carried out around the line Y 1 = 0, where discontinuity occurs (see [16] for more details about the singularity of the K-O problem).We further report the time evolution of the variance of y 1 with = 2, 10 in Figure 7. Little difference is observed for the two simulations.For brevity, we do not report the results for y 2 and y 3 , yet similar observations can be obtained.We remark that the numerical results reported agree well with the benchmarks in the literature [26].

Burgers' Equation with Random Initial Conditions
The last example we consider in this paper is the inviscid one-dimensional Burgers' equation: with the following random initial condition where . Note that there are three states for the initial condition.As time evolves, the solution will develop two shocks with different speed.After some time, the two shocks will interact with each other and finally form one single shock.It is quite challenging to simulate this problem, since this numerical solution is discontinuous in both physical and random space.Below, we will apply the proposed AMR WENO stochastic collocation method to solve this problem.We use a fifth order finite difference WENO scheme coupled with a third order SSP Runge-Kutta method to solve Burgers' equation, which is known as a highly accurate and robust numerical scheme for solving general hyperbolic conservation laws [28].The numerical solution is computed up to T = 0.0543.In Figure 8, we report the contour plots of solution u in the random space with = 1, 5 at location x = 0.79, where the shock interaction occurs.The associated adaptive meshes are also provided.Again, it is observed that the local mesh refinement is carried out around the discontinuities of the solution in the random space.

Conclusions and Discussion
In this paper, we developed an adaptive high-order AMR WENO-based multi-element stochastic collocation method for differential equations with random coefficients and non-smooth solutions, which consists of two important components including a WENO interpolation methodology and an AMR stochastic domain decomposition framework in random space.The numerical experiments were performed for several benchmark tests and the efficacy and efficiency of the proposed method was verified.The scheme is based on a tensor-product construction and hence suffers from the "curse of dimensionality".Future work is to couple the AMR and WENO methodologies in random space with sparse grid [23,26] and adaptive Analysis Of Variance (ANOVA) [11,18] techniques which are efficient tools known for approximating high-dimensional problems.

Figure 7 .
Figure 7.The time evolution of the variance of y 1 .

Table 2 .
Approximation investigation of f 2 .The error of variance (var.err.) for linear reconstruction and WENO reconstruction.