An efficient hybrid method for uncertainty quantification

A technique for coupling an intrusive and non-intrusive uncertainty quantification method is proposed. The intrusive approach uses a combination of polynomial chaos and stochastic Galerkin projection. The non-intrusive method uses numerical integration by combining quadrature rules and the probability density functions of the prescribed uncertainties. A stable coupling procedure between the two methods at an interface is constructed. The efficiency of the hybrid method is exemplified using hyperbolic systems of equations, and verified by numerical experiments.

interfaces is standard practice in many fields of computational mechanics and computational physics. Practical examples include the introduction of mesh-refinement or coarsening to increase accuracy and/or efficiency. By also considering uncertain external data in the boundary conditions and forcing functions, the question arises: can we use different uncertainty quantification (UQ) methods in different regions of the computational domain separated by interfaces?
In [18] it was shown that intrusive polynomial chaos with stochastic Galerkin projection (PC) is more efficient than non-intrusive numerical integration (NI) for slowly varying (in stochastic space) stochastic problems, while NI is more efficient for rapidly varying problems. This result suggests that a combination of the two techniques could be efficient. These methods must of course be coupled in an accurate and stable way, and this is the topic of this paper.
We aim for a method which couples an intrusive and non-intrusive UQ method. The intrusive method combines a polynomial chaos basis with stochastic Galerkin projection [4,10,19]. The non-intrusive method uses the combination of a quadrature rule and the probability density functions of the uncertainties prescribed [1,3,15]. The proposed interface coupling for the different UQ methods is proven to be stable, by extending the work in [5], for non-conforming spatial finite difference methods [2,16,17]. Our intention is to exploit the different properties of intrusive and nonintrusive UQ methods and obtain a higher efficiency by combining them.
The rest of the paper proceeds as follows. The hyperbolic system of equations is introduced in Sect. 2. Section 3 presents the intrusive method and Sect. 4, the non-intrusive one. Section 5 describes the interface coupling treatment. A stable and accurate finite difference formulation is presented in Sect. 6. The theoretical results are confirmed using numerical experiments in Sect. 7. Finally, conclusions are drawn in Sect. 8.

The continuous problem
To present the technique it suffices to consider a hyperbolic system of equations posed on two adjacent spatial domains (Ω L and Ω R ) x ∈Ω L , t > 0, (2.1a) x ∈Ω R , t > 0, (2.1d) The solutions on the left and right domains are denoted by u = [u 1 (x, t, ξ), . . . , u M (x, t, ξ)] and v = [v 1 (x, t, ξ), . . . , v M (x, t, ξ)], respectively. Moreover, the random variable ξ models the uncertainty. The boundary operators defined on the boundaries ∂Ω L and ∂Ω R are denoted L L and L R , while f L (x, ξ), f R (x, ξ) and g L (x, t, ξ), g R (x, t, ξ) are the stochastic initial and boundary data, respectively. The data is assumed to be fully compatible at the boundaries and interface. Further, we let Ω L , Ω R and Γ represent the left and right domains, and interface boundary, respectively, see Fig. 1. Further, we consider the M × M matrix A to be where + A and − A are positive and negative definite diagonal matrices (In the case where the matrix A is not diagonal, we assume that it can be symmetrized and diagonalized).

Well-posedness
Applying the energy method to (2.1) (that is multiply (2.1a) and (2.1d) by u and v, respectively, and integrate over the particular spatial domain and adding them together), considering only the terms at the interface (x = 0) and using the interface condition u x=0 = v x=0 , yields We can conclude that the energy in (2.3) is preserved and that proper outer boundary conditions leads to well-posedness.

Polynomial chaos
In the polynomial chaos method, the solution is expanded in a polynomial basis, In (3.1), u k denotes a vector of expansion coefficients. We choose the polynomial basis functions {ψ(ξ )} ∞ i=0 to be orthonormal with respect to the inner product In (3.2), P(ξ ) is the probability measure of ξ defined on the stochastic domain Ω ξ . The coefficients u k in (3.1) are given by the projection The mean can be expressed in terms of the coefficients as since E[ψ 0 ] = 1 and E[ψ k ] = 0, ∀k > 0. Moreover, the variance is given by Further details of the PC framework can be found in [14].

Stochastic Galerkin projection
In order to apply the stochastic Galerkin projection method in the left domain Ω L , we start by truncating the series (3.1) as where a slight abuse of notation has been used. Inserting (3.6) into the problem in the where v is the solution from the right domain (Ω R ). By multiplying (3.7) with ψ l , for l = 0, 1, . . . , M PC and integrating over the stochastic domain Ω ξ , we obtain

Numerical integration
Numerical integration in one dimension is formulated as where f is the function being integrated and ρ is the density function. The integer M N I denotes the number of quadrature points ξ m , and w m denotes the quadrature weights. Extending (4.1) to higher dimensions can be done in a straightforward manner as where Ω ξ is the p-dimensional domain andξ = (ξ 1 , . . . , ξ p ). The quadrature points and weights for dimension i are denoted ξ m i and w m i respectively. Adaptive sparse grid techniques for high-dimensional problems [1] can be used to improve efficiency and ease the "curse of dimensionality" of the quadrature rule. In this work however, the 4th-order accurate Simpson's rule [3] is used for all computations.

Interface treatment
To couple the problem (2.1) in the stochastic dimension at an interface in space, a stable interface treatment is required. The specific interface treatment we consider applies and extends the work done in [5]. In this technique, an intermediate so-called glue grid G h (with a corresponding glue space G ), together with a set of projection operators is introduced. The solutions at the interface are projected onto the glue grid where their difference is penalized, leading to a provably stable method. The glue space consists of piecewise polynomials of a sufficiently high degree. We choose Legendre polynomials since they form an orthogonal basis with respect to the uniform density function on the interval (−1, 1). Further, we letF l2g andF g2l be the projection operators between the left grid l and the glue grid G and vice verse, respectively. The operators projecting values on the right grid to the glue grid and back are denotedF r 2g andF g2r . The projection operators are constructed to satisfy the following inner product preserving relation [5,11] where M = γ ψ(ξ )ψ T (ξ ) dξ is a symmetric positive semi-definite matrix and the ψ's are basis functions of the glue grid. The matrices P ξ L and P ξ R represent suitable quadrature rules on the left and right side of the interface. Further, the projection operators are required to satisfy a set of accuracy conditions, namely that the glue grid can represent qth order polynomials exactly. This means that a polynomial (in ξ ) of order q projected onto a grid is mapped exactly onto the glue grid byF l2g (orF r 2g ). The order q is chosen based on the order of accuracy of the uncertainty quantification methods used on the left and right domain, see below for details.

Remark 5.1
Note that we do not requireF l2gFg2l = I orF r 2gFg2r = I . This means that a projection from the left (or right) to the glue grid and back is not required to be exact.
The reader is referred to [5] for a detailed description of the projection operators. Next, we briefly describe the cases that we will consider.
-Coupling: NI to NI The glue grid and projection operators are analogous to the ones in [5]. A brief presentation of these operators and the corresponding glue grid technique can be found in Appendix A.1.

-Coupling: PC to PC
When coupling PC to PC, the glue grid consists of a polynomial basis of the same degree as the maximum degree on the left and right domain (in order to satisfy the accuracy conditions). The projection operators are constructed in a straightforward manner, projecting a set of polynomials to another set (see Appendix A.2 for a simple example). -Coupling: NI to PC The coupling between NI and PC is done in a similar way as in the coupling between NI and NI. The number of polynomial basis functions per subinterval in the glue grid is chosen to be same as the number of basis functions of the PC grid to maintain the accuracy conditions in the coupling. For more details, see Appendix 3.

The semi-discrete formulation
The numerical approximation of (2.1) based on summation-by-parts (SBP) operators with simultaneous approximation terms (SAT) [2,8,12,17] is where To impose the boundary and interface conditions the so-called SAT technique is used. The terms that contains the 's are the so-called SAT penalty terms, which forces the solution to satisfy the boundary and interface conditions weakly. The penalty matrices ( 0 ) L and ( N ) R will be chosen such that a stable spatial boundary treatment in the left and right domain is obtained. Furthermore, ( N ) L and ( 0 ) R are used to enforce continuity of the solution at the interface. Finally, ( I ) L,R are correction penalty matrices used to cancel the interface contributions in stochastic space coming from the operators D L,R . To clarify: the indices (0, I , N ) represent the first, the interface, and the last grid points in space, respectively. Similarly, L and R represent the right and left domain.
The numerical solution on the left sub-domain U is arranged as and the numerical solution V in the right sub-domain is organized in the same way.
where m denotes the vector component and j the PC coefficient. The data is organized depending on the stochastic method as . . .
for PC and for NI. In (6.3) and (6.4),f andḡ are the initial and boundary vectors as a function of ξ , respectively.

Stability
The discrete energy method is applied to (6.1) (we multiply the equations in (6.1) by U T (P L ⊗ P ξ L ⊗ I M ) and V T (P R ⊗ P ξ R ⊗ I M ), respectively, add them together, and use the SBP property (Q L,R are almost skew-symmetric)). This yields when letting P ξ L,R = (P ξ L,R ⊗ I M ) and ignoring the outer boundary terms as in Sect. 2.1. The norms used in (6.5) are U 2 Further, using (5.1) in (6.5), (for example the term 2U T where and we have decomposed the SAT penalty terms as In (6.6), the first two terms on the right-hand side come from the correction terms (the difference between U or V and its projection to the glue grid and back) in (6.1), while the third term contains the interface terms on the glue grid.
To cancel the interface terms not involving the operatorsF ·2· we let (¯ I ) L = −A/2 and (¯ I ) R = A/2, (other choices are also possible, see for example [13]) which results in d dt which implies that the energy is preserved (ignoring the outer boundary treatment) just as in (2.3).

Remark 6.1
The quadratures P ξ L and P ξ R can be replaced by any other positive definite matrix without loss of stability, however, to ensure accuracy, the operators F ·2g and F g2· have to satisfy certain accuracy conditions. These accuracy conditions involve projecting polynomials of a specific degree. In all our computations, we have, as in [5] used SBP quadratures of the same accuracy as the number of polynomial basis functions of the glue space.

Numerical experiments
In this section we will illustrate the coupling technique by considering applications in one and two dimensions.

Advection in one dimension
We consider (2.1) with As a measure of comparison, the quantity is used. In (7.2), U is the computed numerical solution (on both domains) and U a a high resolution (in ξ ) numerical solution on the same spatial deterministic grid as U . In the computations below, 5th order SBP-operators with 50 grid points in each domain and the classical 4th order Runge-Kutta as time integrator are used when computing U and U a . The uncertainty ξ is uniformly distributed between −1 and 1. Finally, the spatial domains we consider are Ω L = (−1, 0) and Ω R = (0, 1), hence the interface is located at x = 0.

Accuracy and efficiency
The order of accuracy p N I when computing the variance (7.2) numerically, for the case of coupling NI with NI is measured as The variance in Var U X ,Y is computed using NI with X and Y number of stochastic grid points in the left domain and right domain, respectively. The constant m denotes the refinement factor of the stochastic grids. To reduce the influence of the deterministic errors, a high resolution solution using 200 stochastic grid points on a single deterministic domain using the same number of deterministic grid points was used. Initial and boundary data and a forcing function was extracted from the following manufactured solution The order of accuracy is computed and shown in the one dimensional case, results in multiple dimensions are analogous. As can be seen from Table 1, the order of accuracy converges to 4, in line with the accuracy of the 4th order Simpson's rule. The spatial operators and time integrator are 5th and 4th order accurate, respectively. The gains in efficiency and accuracy that we will show computationally below, rely on the cost of the glue grid procedure being cheap. The glue grid is only applied at the interface, hence contributes linearly to the total computational cost (which overall is quadratic) with respect to the number of stochastic quadrature points and PC coefficients. The computational cost of the interface treatment is hence small compared to the overall computational cost.

Numerical restults: coupling NI to NI
To exemplify the potential efficiency gain in the coupling procedure, the problem (2.1) using a manufactured solution is constructed in such a way that the uncertainty is varying slowly in the left domain compared to the right. We choose The manufactured solution (7.5) contains the term 1/(1 + e −0.1ξ ) which is a S-shaped function ranging from 0 to 1 or as in this particular case between 0 and the factor e −100(x−1/2) 2 . A sine function is then added to produce variations. A smaller factor (as in the left domain (−1 < x < 0)) gives a slower varying solution and larger factor gives a faster varying solution (as in the right domain (0 < x < 1)). As before we use 200 stochastic grid points for our reference solution. Figure 2 shows the error of the variance as a function of the total number of stochastic grid points. As can be seen, a higher number of grid points in the right (more rapidly varying) domain reduces the error in the variance compared to the other cases. Since all cases have the same total number of grid points, we can conclude that having a higher number of grid points in the right (more rapidly varying) part of the domain is also more efficient.

Numerical results: coupling PC to PC
The coupling between PC and PC is exemplified using the manufactured solution (7.5). As a reference, a solution using 25 PC coefficients in each domain is used. Figure 3 illustrates the error as a function of the total number of basis functions. Similarly to the NI-NI case above, a higher number of PC coefficients in the right (more rapidly varying) domain results in a lower error in the variance. Similarly to the case above, we can conclude that having a higher number of PC coefficients in the right (more rapidly varying) part of the domain is the most efficient.

Numerical results: coupling NI to PC
To illustrate the coupling between NI and PC, the manufactured solution (7.5) is used. Figure 4 shows the error in the variance as a function of the number of stochastic grid points (for NI) and basis functions (for PC). In the comparison, we increase the number of stochastic grid points M N I with the number of basis functions as M N I (M PC ) = 30(M PC − 3), e.g. if M PC = 5, then M N I = 60. As can be seen, the combination of using NI on the rapidly varying part (in terms of ξ ) and PC on the slowly varying part of the domain performs better than the reverse combination. This is an indication that combining NI and PC could outperform the use of the same UQ method in the whole domain. The results in Fig. 4 shows an error reduction by approximately one order of magnitude, which show that choosing a suitable UQ method in certain regions is important. The most efficient choice is having NI in the rapidly varying part and PC in the slowly varying part of the domain.

Maxwell's equations
To exemplify the technique in multiple dimensions, we consider Maxwell's equations in two dimensions [17]. The equations relating the electric and magnetic fields are and μ correspond to the electric field, magnetic field, electric current density and charge density, permittivity and permeability. We simplify (7.6) by letting J = 0 and ρ = 0 and rewrite it on matrix form as Su t + Au x + Bu y = 0. (7.7) The matrices in (7.7) are with the solution u = H z , E x , E y T . Similar to the one dimensional case, we have an interface at x = 0, 0 ≤ y ≤ 1. For simplicity, the boundary conditions are chosen to be the characteristic ones, i.e. we specify In the numerical experiments we let = 1 and μ = 1 be defined on the spatial domains Ω L = (−1, 0) × (0, 1) and Ω R = (0, 1) × (0, 1). The coupling is exemplified using the manufactured solution which varies more stochastically in the left domain compared to the right.

Coupling NI to NI
Numerical results when coupling NI with NI were computed using a reference solution calculated using 150 stochastic grid points in the left and right domain using NI in both domains. The reference solution is computed using an interface in order to reduce the deterministic errors in the comparison. Figure 5 illustrates the error in the variance with respect to the total number of stochastic grid points. As in the one dimensional case, we note that a higher number of grid points in the more rapidly varying (left) part of the domain gives the smallest total error in the variance, as well as a higher efficiency.

Coupling PC to PC
Coupling PC with PC was exemplified using a reference solution computed with 25 PC coefficients in both domains. The numerical results are depicted in Fig. 6, where the  The error in the variance as a function of the number of PC basis functions when coupling NI with PC error in the variance is shown as a function of the number of PC coefficients. Again, we note that a higher number of points in the more rapidly varying (left) domain reduces the error in the variance the most. As before we can conclude that the above mentioned case is the most efficient one.

Coupling NI to PC
The numerical results when coupling NI with PC were performed with a reference solution computed using 150 stochastic grid points in the left and right domain using NI in both domains. Figure 7 shows the error in the variance as a function of stochastic grid points (in NI) and PC coefficients (in PC). Results when combining PC on the left with NI on the right domain and vice verse are illustrated. We note that, as in the one dimensional case, having NI in the rapidly varying domain and PC in the slowly varying domain outperforms the reverse combination. The reduction of errors in the variance is approximately two orders of magnitude. Finally, we can conclude that the most efficient case coincides with the results in the one dimensional setting, i.e. having NI in the rapidly varying domain and PC in the slowly varying domain.

Conclusions
A stable, accurate interface coupling between the uncertainty quantification methods NI-NI, PC-PC and NI-PC has been constructed. The coupling procedure consisted of constructing an intermediate so-called glue grid, onto where the solutions at the interface were projected and penalized using corresponding projection operators. The coupling procedure was exemplified using hyperbolic systems of equations in one and two dimensions.
Numerical results when coupling NI-NI and PC-PC clearly showed that a higher proportion of grid points/coefficients in a more rapidly varying part of the domain gives the smallest error in the variance. Coupling NI-PC gave better results when having NI in the more rapidly varying domain, resulting in a reduction of the errors in the variance of 1-2 orders of magnitude compared to the reverse combination. The performance benefits shown by the numerical experiments when using different UQ methods in different regions, gives a strong indication that combining different UQ techniques can be efficient for certain types of problems.
To extend the technique to multiple stochastic dimensions, is in principle straightforward, but some technical problems requires further research. One has to construct stable interface couplings in multiple dimensions, which can be done following the guidelines in [6,7,9].
Funding Open access funding provided by Linköping University.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A projection operators
This section gives a more detailed description of the various projection operators.

A.1 Description of projection operators coupling NI with NI
The order of the polynomial basis (number of basis functions per subinterval) of Legendre polynomials (i.e. the basis of the glue grid) is set to equal the accuracy order of the NI method (denoted q N I ) in order to satisfy the accuracy condition. As an example, Fig. 8 shows a simple illustration when coupling NI with NI using six and four stochastic grid points on the left and right domain, respectively. The intervals of the glue grid is depicted in Fig. 8. The projection operators from the left grid to the glue grid and vice versa in this example (using q N I = 2), are (A.1) Projection operators from the right grid to glue grid have similar structures as the ones in (A.1), hence will not be shown here.

A.2 Example of projection operators coupling PC with PC
As an example, we consider a coupling between four (on the left domain) respectively two (on the right domain) basis functions, see Fig. 9. The glue space is in this case a space of polynomials of degree 4 − 1 = 3 (the maximum degree of the PC expansions on the left and right domain). The projection operators are in this example The coupling between PC and PC with four and two basis functions respectively with an intermediate glue grid Note that the operators in (A.2) satisfy F g2l F l2g = I 4 and F g2r F r 2g = I 2 , where I 2,4 are identity matrices of dimension two and four respectively.

A.3 Description of projection operators coupling NI with PC
The number of polynomial basis functions per subinterval in the glue space is chosen to be same as the number of PC coefficients (or the order of the NI method (q N I ) if it is greater), in order to maintain the accuracy conditions of the coupling. The projection operators from the NI grid to the glue grid are hence identical to the ones described in the section Coupling: NI to NI. The projection operators between PC and the glue grid and back are transformations between a set of basis functions (polynomials) to another on different intervals (the intervals are determined by the NI grid). This means that coefficients corresponding to basis functions on the whole interval are simply mapped to coefficients corresponding to basis functions on the different intervals (which constitutes the glue grid). As an example of the transformation between the PC grid and glue grid, Fig. 11 illustrates a basis function on the PC grid and Fig. 12 the corresponding function projected on the different intervals in the glue grid. The polynomial basis used in this work are the Legendre polynomials due to the uniform distribution of the uncertainty. As an illustrative example, Fig. 10 shows the coupling between NI and PC using five and four stochastic grid points/basis functions on the left (NI) and right (PC), respectively.