A DYNAMIC BI-ORTHOGONALITY BASED APPROACH FOR UNCERTAINTY QUANTIFICATION OF STOCHASTIC SYSTEMS WITH DISCONTINUITIESThis work was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2010-0025484)

A Dynamic Bi-Orthogonality Based Approach for Uncertainty Quantification of Stochastic
Systems With Discontinuitiesthanks: This work was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2010-0025484)

Piyush M. Tagade Post-doctoral Research Fellow, Division of Aerospace Engineering, Korea Advanced Institute of Science and Technology, Republic of Korea. (piyush.tagade@kaist.ac.kr).    Han-Lim Choi Assistant Professor, Division of Aerospace Engineering, Korea Advanced Institute of Science and Technology, Republic of Korea. (hanlimc@kaist.ac.kr). Questions, comments, or corrections to this document may be directed to that email address.

The use of spectral projection based methods for simulation of a stochastic system with discontinuous solution exhibits the Gibbs phenomenon, which is characterized by oscillations near discontinuities. This paper investigates a dynamic bi-orthogonality based approach with appropriate post-processing for mitigating the effects of the Gibbs phenomenon. The proposed approach uses spectral decomposition of the spatial and stochastic fields in appropriate orthogonal bases, while the dynamic orthogonality condition is used to derive the resultant closed form evolution equations. The orthogonal decomposition of the spatial field is exploited to propose a Gegenbauer reprojection based post-processing approach, where the orthogonal bases in spatial dimension are reprojected on the Gegenbauer polynomials in the domain of analyticity. The resultant spectral expansion in Gegenbauer series is shown to mitigate the Gibbs phenomenon. Efficacy of the proposed method is demonstrated for simulation of a one-dimensional stochastic Burgers equation with uncertain initial condition.

Key words.    Uncertainty Quantification, Gibbs Phenomenon, Stochastic Spectral Methods, Dynamically Bi-orthogonal Field Equations

AMS subject classifications.    74S25, 74S60, 74J40

1 Introduction

Continuous advancements in digital technologies have resulted in ubiquitous use of computer simulators for investigation of many large-scale systems. The simulator predictions are often uncertain due to unknown/poorly known physics, model parameters, initial and boundary conditions. Need and importance of uncertainty quantification in simulation predictions has already been emphasized by researchers in varied fields [34, 28, 29, 33]. Subsequently, research community have invested significant efforts to develop various uncertainty quantification methodologies (see [45] and references therein.). In essence, uncertainty quantification is an estimate of joint probability distribution of system responses conditional on probabilistic specification of uncertainties in the simulator model, parameters, initial and boundary conditions etc.

1.1 Background

For further discussion, consider a probability space where is a set of elementary events, is associated -algebra and is a probability measure defined over . A generic stochastic partial differential equation (SPDE) on is given by


where is a random variable, is time, is a spatial dimension and is an arbitrary differential operator of the form


being flux. The initial condition of the system is specified using a random field , while, the boundary conditions are given by


where is a linear differential operator. Without loss of generality, the discussion in this paper is presented assuming to be a scalar function defined over a one-dimensional space, however, the discussion can be extended to a more generic case in straightforward manner.

Owing to generality and ease of implementation, Monte Carlo methods are widely used for solution of (LABEL:spde) [37]. However, the Monte Carlo method typically requires a large number of samples for acceptable accuracy [37, 5]. For simulation of large scale systems, computational cost may render application of Monte Carlo methods for uncertainty quantification intractable. Stochastic spectral projection (SSP) based methods can provide a computationally efficient alternative to Monte Carlo methods with comparable accuracy. Under a generic condition of square integrability, can be represented as a parametric field in a Hilbert space  [1]. The SSP methods use the tensor product representation of in appropriate Hilbert spaces and with possible choices  [44]


A Karhunnen-Loeve expansion [12] uses the tensor product (LABEL:tens_prod_1), where is spectrally represented in terms of orthogonal basis in with temporally varying expansion coefficients in . A generalized Polynomial Chaos (gPC) expansion [50] is an example of tensor product (LABEL:tens_prod_3) with orthogonal basis in and the associated expansion coefficients in . This paper investigates a bi-orthogonal method, which uses the gPC basis in , while, dynamically orthogonal eigenfunction basis is used in .

One of the earliest exposition of a SSP-based method is the Homogeneous Chaos theory introduced by Wiener [47, 48], where random variables are expanded in terms of orthogonal Hermite polynomials. The expansion converges in mean square sense for any nonlinear functional [2]. Meecham and Jeng [27] used Homogeneous Chaos for turbulent modeling, however it was found that the convergence of chaos expansion is slow for turbulent field [38, 6]. Ghanem and Spanos [11, 10, 12] have proposed Polynomial Chaos (PC) method based on Homogeneous Chaos theory for uncertainty propagation in structural mechanics simulations. Subsequently, the method has been widely used for uncertainty propagation in fluid flow simulations [21, 24, 26, 32]. Present state of the art for SSP methods is based on the Galerkin projection based generalized Polynomial Chaos method introduced by Xiu and Karniadakis [49, 50], which uses the representation of form (LABEL:tens_prod_1). Using separability, is decomposed in orthogonal polynomial chaos basis , where belong to the Askey scheme of orthogonal polynomials [19] with weight function given by probability distribution function of an uncertain input. Thus, using gPC, the dynamical response is spectrally represented as


where are dynamical expansion coefficients. See Mathelin et al. [25] and Najm [30] for extensive review of the research work in this field.

Remark 1

Throughout this paper, the operator is used interchangeably to represent ‘equal to’ and ‘approximately equal to’. Note that all the instances of in spectral representation defines the ‘approximately equal to’ operator.

1.2 Motivation

Research work presented in this paper is motivated by the need to quantify uncertainty in simulation of systems with discontinuous solutions. In particular, the paper deals with stochastic systems defined using the hyperbolic conservation laws that invariably results in discontinuous solution irrespective of the initial conditions. Simulation of such system have attracted significant interest by research community and literature is rich with methods for resolution of discontinuities in deterministic simulations (see [22] and references therein for review of the methods).

The use of spectral methods to hyperbolic SPDEs is known to pose following significant numerical challanges [31]

  1. numerical solution using finite difference schemes lead to spurious oscillations in expansion coefficients and/or orthogonal basis

  2. spectral expansion of the solution of hyperbolic SPDEs lead to Gibbs phenomenon as discontinuities develop in solution [31]. The Gibbs phenomenon is characterized by [17, 14] a) slow convergence of spectral approximation of the solution at points away from the discontinuity; and b) oscillations are observed in the solution near discontinuities that do not decrease with increasing .

There are recent research efforts that focus on addressing the issue of Gibbs phenomenon in gPC settings [31]. Wan and Karniadakis [46] have proposed a multi-element gPC (ME-gPC) method to mitigate the effect of Gibbs phenomenon (also see Lin et al. [23]). The ME-gPC method uses domain decomposition of the stochastic space with locally defined gPC basis, while the forward problem is solved for each subdomain. However, the computational cost of the ME-gPC method rises quickly as the number of subdomains increases. Poette et al. [36] have proposed an entropy based intrusive polynomial moment method for uncertainty propagation in non-linear systems with discontinuous solutions. The method reformulates the SPDE in terms of appropriately selected entropic variables that are bijection of the uncertain parameters, while, the gPC expansion of the entropic variables is defined through the Galerkin projection of the bijection. The method is implemented in two steps: in the first step, the resultant stochastic Galerkin system is discretized using a finite volume approach, and upwinded Row scheme is used for numerical solution. In the second step, gPC expansion of the entropic variables is calculated so that the entropy of the system is minimized. The method is demonstrated for solution of the stochastic Burgers equation. The method is computationally efficient than the ME-gPC, however, the requirement for minimization of the entropy at each time step results in the computational cost higher than the gPC method. Tryoen et al. [43] uses Roe-type solver for intrusive Galerkin projection of uncertain hyperbolic systems. The solver uses Dubois and Mehlman entropy corrector that avoid certain entropy violating solutions. Sargsyan et al. [40] have proposed a two step method for uncertainty quantification of systems with discontinuous response. In the first step, a limited number of simulation runs are used to infer the shock location using the Bayesian inference. In the second step, localized gPC basis is defined in the region of analiticity, while intrusive Galerkin projection is used to propagate the uncertainty. In an alternate non gPC setting, Chantrasmi et al. [4] have proposed a Pade-Legendre interpolant based approach, where, simulation runs are obtained at predefined Gauss-Legendre quadrature nodes and the resultant discontinuous system response is reconstructed using the Pade interpolation. Despite of these recent progresses, resolution of Gibbs phenomenon in application of SSP methods to stochastic hyperbolic systems is an area of current research interest [40].

1.3 Proposed Method

Existence and resolution of Gibbs phenomenon for deterministic functions is extensively investigated in the literature [17, 14, 16]. Spectral expansion of a function using global orthogonal basis is contaminated by the presence of discontinuities. However, expansion coefficients contain enough information to recover the function with high accuracy using appropriate post-processing. Gibbs phenomenon can be completely resolved by reprojecting the partial sum on Gibbs complementary basis [14]. Gottlieb et al. [15, 13, 14, 16] have shown that exponentially convergent approximation can be obtained at point values in subinterval of analyticity by reprojecting the partial sum on orthogonal Gegenbauer polynomial basis. The post-processing method can be extended to stochastic functions, provided, the spectral expansion is available in terms of orthogonal basis in spatial dimensions.

Note that the gPC expansion (LABEL:gpc) uses orthogonal basis in while corresponding coefficients are functions in , though, the discontinuity reside in spatial dimension. Thus, enough information is not available in to resolve Gibbs phenomenon using Gegenbaeur reprojection method. Resolution of Gibbs phenomenon proposed in this paper is based on the observation that exponential accuracy can be recovered if is expanded in terms of orthogonal bases in spatial dimension, and subsequently reprojected on an approximate Gibbs complementary basis.

Tagade and Choi [42, 41] have proposed a dynamic bi-orthogonal field equations (DBFE) method for solution of SPDEs in the context of Bayesian calibration. This paper proposes a method for solution of hyperbolic SPDEs using DBFE, that exploits orthogonal decomposition of the spatial dimension to develop a post-processing approach for mitigation of the Gibbs phenomenon. Consider a generic Karhunnen-Loeve expansion of


where is the mean, are eigenfunctions which form complete orthogonal basis in at a given time step , while are independent zero-mean random variables. Using polynomial chaos expansion of in (LABEL:kl_expn), bi-orthogonal expansion of is given by


A Dynamic Orthogonality (DO) condition [39] is used to derive closed form solution of evolution equations for , and . However, the resultant approximation of in (LABEL:kl_expn_dbfe) using the solution of , and exhibits the Gibbs phenomenon. A Gegenbauer reprojection based post-processing method is proposed to mitigate the effect of Gibbs phenomenon by reprojecting the eigenfunctions on the Gegenbauer basis. The proposed method is demonstrated for solution of a stochastic one-dimensional Burgers equation with uncertain initial conditions.

The rest of the paper is organized as follows. In section 2, mathematical formulation of the method is provided. Proposed post-processing method for Resolution of the Gibbs phenomenon using dynamic bi-orthogonality based method is discussed in section 3. Section 4 provides numerical results for Burgers equation to demonstrate efficacy of the proposed method. Finally, paper is summarized and concluded in section 5.

2 Dynamically Bi-orthogonal Field Equations (DBFE)

Definition of the following operators will be used in the remainder of the paper.

Definition 1

Inner product on for a given is defined as


Similarly, inner product on is defined as


The expectation operator on is defined as


Further using definition of the expectation operator, the covariance is defined as


The DBFE method uses the bi-orthogonal expansion (LABEL:kl_expn_dbfe) of in (LABEL:spde). However, independent evolution equations for the mean , the eigenfunctions and the corresponding coefficients can not be obtained concurrently using the SPDE in (LABEL:spde). Well posed evolution equations for the unknown quantities can de derived by imposing the additional dynamic orthogonality (DO) condition [39]. The DO condition is specified as [39]


where is the number of eigenfunctions retained in the expansion (LABEL:kl_expn_dbfe). Note that the DO condition constrains the dynamic evolution of the eigenfunction basis such that the orthonormality of the eigenfield is retained at all time steps [39].

Using the DO condition and bi-orthogonal expansion of , SPDE (LABEL:spde) can be reformulated into a set of PDEs and ODEs as follows.

Proposition 1

Using the DO condition, dynamic evolution equations of , and are given by


Proof. Proof of (LABEL:mean_ev) and (LABEL:efun_ev) is given by Sapsis and Lermusiaux [39], which is briefly presented here for completeness, while, (LABEL:eyp_ev) is derived here by introducing the bi-orthogonality.

Proof of (LABEL:mean_ev) and (LABEL:efun_ev)

Use a generic KL expansion (LABEL:kl_expn) in (LABEL:spde) to obtain


Application of the expectation operator to (LABEL:spde_kl) gives (LABEL:mean_ev).

Multiply (LABEL:spde_kl) by and apply the expectation operator to obtain


Multiply (LABEL:cov_dyidt) by , take inner product and apply the DO condition to get


Use (LABEL:cov_dyidtuk) in (LABEL:cov_dyidt) to obtain


which can be written in the matrix as


where is the covariance matrix with element .

Multiply both sides of (LABEL:spde_kl) by , take inner product and use the DO condition and orthonormality of eigenfunctions to get


which on application of expectation operator gives


Using (LABEL:inn_prod_uij_2) in (LABEL:inn_prod_uij_1) gives

Proof of (LABEL:eyp_ev)

Basic conditions of KL expansion ensures that are square integrable random variables, thus, according to Cameron and Martin theorem [2], can be approximated to arbitrary accuracy using spectral expansion in terms of orthogonal basis in . Hence, can be spectrally represented using terms as [50]


where are orthogonal basis in . Differentiating (LABEL:pc_exp) with respect to gives


which on Galerkin projection provide


Having (LABEL:bi_eq2_1) and (LABEL:bi_eq2_3) in (LABEL:gov_eq2) results in (LABEL:eyp_ev).


2.1 Initial and Boundary Conditions for DBFE

The initial conditions for DBFE are defined by specifying the mean , the eigenfunctions and the coefficients , which are given by considering the bi-orthogonal expansion of the initial condition of SPDE (LABEL:spde) as


Applying the expectation operator, initial condition for the mean field is given by


The initial condition for eigenfunctions, , is given by solution of the Fredholms’ integral equation of the second kind [3]


where is the covariance function of and are the eigenvalues. A Galerkin projection based method is used to solve (LABEL:fie) numerically [18]. For a Gaussian process, all the coefficients are zero except


thus defining as normal variables with variance . For a generic non-Gaussian case, the initial expansion coefficients are given by


Boundary conditions for DBFE are given by using the generic KL expansion of , which specifies the boundary conditions for the SPDE in (LABEL:spde), as


Applying the expectation operator on (LABEL:bound_kl_bi), boundary condition for the mean field is given by


Multiply (LABEL:bound_mean) by and apply the expectation operator


which can be specified in a matrix form as


where is the covariance matrix.

3 Resolution of the Gibbs Phenomena using Dynamically Bi-orthogonal Field Equations

Key challenges involved in the numerical implementation of the DBFE for SPDEs with discontinuous solutions are: (a) numerical evaluation of the mean and the eigenfunctions (using (LABEL:mean_ev) and (LABEL:efun_ev)) results in the development of spurious oscillations as discontinuities evolve in the solution; (b) the resultant bi-orthogonal expansion (LABEL:kl_expn_dbfe) exhibits the Gibbs phenomenon characterized by the oscillations near the discontinuity location and the slow convergence away from the discontinuity location. This paper proposes a two step approach for the application of DBFE to the stochastic systems with discontinuous solution, that exploits the orthogonal decomposition of the spatial field in the DBFE approach. In the first step, extension of the existing schemes is proposed to derive the non-oscillatory numerical scheme for DBFE. In the second step, a Gegenbauer reprojection based method is proposed to mitigate the effects of the Gibbs phenomenon. In this section, the proposed approach is described in detail.

3.1 Numerical Scheme for Non-oscillatory Solution

Due to discontinuities, numerical solution of (LABEL:mean_ev)–(LABEL:eyp_ev) is expected to contain spurious oscillations with reduced accuracy in and . Total Variation Diminishing (TVD) finite difference schemes are widely used in the literature to obtain non-oscillatory accurate weak solutions to hyperbolic partial differential equations [22]. This paper demonstrates extension of a TVD scheme for numerical solution of (LABEL:mean_ev)–(LABEL:eyp_ev), however, note that other numerical schemes can similarly be extended for DBFE.

For notational convenience, this subsection defines as an approximate value of at the grid point and present time . Consider a numerical scheme applied to SPDE (LABEL:spde) for a given as


where is an interpolated flux inside a cell . A point numerical scheme is defined by using points on either side of for interpolation, thus,


such that [22]


is also defined in a similar manner. The TVD scheme (LABEL:tvd_1) can be extended to DBFE by using


for numerical solution of (LABEL:mean_ev)–(LABEL:eyp_ev). The resultant numerical scheme retain non-oscillatory property in and .

In the present paper, the DBFE method is implemented using the first order central scheme proposed by Kurganov and Tadmor [20] (central KT scheme), for which numerical flux is defined as


where is the maximum local speed at . The gPC expansion of is given by


where are gPC expansion coefficients. Use (LABEL:gpc_spec_rad) in numerical flux (LABEL:num_flux) to obtain


where and .

The numerical flux (LABEL:num_flux_dbfe) is used in (LABEL:tvd_1) to obtain the first order central difference scheme for numerical solution of (LABEL:mean_ev)-(LABEL:eyp_ev). order Runge-Kutta method is used for time integration in (LABEL:mean_ev) and (LABEL:efun_ev), while Euler method is used to solve ODE (LABEL:eyp_ev). Gauss-Legendre quadrature is used to calculate inner products in spatial dimensions, while, Monte Carlo integration is used for inner products in stochastic dimensions.

3.2 Post-processing using Gegenbauer Polynomial

Let the solution of the SPDE in (LABEL:spde) is obtained using the DBFE method till time such that the discontinuities are developed. For a given , let denote an interval of analyticity for , where a random variable can be defined such that . Gibbs phenomenon can be resolved in the interval by reprojecting the partial sum (LABEL:kl_expn_dbfe) on a suitable Gibbs complementary polynomial basis [16]. In the present paper, the partial sum is reprojected on a Gegenbauer Polynomial, which is defined for a as follows.

Definition 2

[16] The Gegenbauer polynomial is a polynomial of degree which is orthogonal over with weight function for . The orthogonality is given by






can be estimated using a recurrence relationship


Using orthogonality of Gegenbauer polynomial, can be represented in terms of exponentially convergent Gegenbauer expansion series truncated at terms as


where ; and . Coefficients are given by


However, is not available for calculation of . Using bi-orthogonal expansion (LABEL:kl_expn_dbfe) in (LABEL:gegen_coef_u),  approximate Gegenbauer expansion coefficients are given by


Note that (LABEL:gegen_coef_g) uses a bi-orthogonal expansion of , thus, is an approximation of . can be expanded in gPC basis as




Coefficients are calculated using Gaussian quadrature by evaluating (LABEL:gegen_coef_g) at collocation points, while, the integrals in (LABEL:gegen_coef_g) are calculated using Gauss-Gegenbauer quadrature. Using , is approximated as


The partial sum (LABEL:gegen_part_sum_pc) approximates with exponential accuracy in the subinterval . Note that determination of domain requires location of discontinuity. An edge detection method [8, 9] can be used to locate points of discontinuities.

4 Numerical Example: 1D Burgers Equation

Efficacy of the proposed method is investigated for solution of a stochastic one-dimensional Burgers equation with uncertain input conditions. Note that previous work has utilized deterministic Burgers equation extensively [51], while, there have been recent investigation on solution of stochastic Burgers equation [35, 36]. Coupled with the ease of implementation, stochastic Burgers equation provide an appropriate context to investigate the efficacy of the proposed method for mitigation of the Gibbs phenomenon.

A one dimensional inviscid stochastic Burgers equation is given by


where . In the present paper, the stochastic Burgers equation is investigated for uncertain initial condition given by


where is a Gaussian process with mean


where and are user defined constants, and the covariance function


where is standard deviation and is the correlation length. In this paper results are presented for . The boundary condition is specified using


4.1 Monte Carlo Method for Stochastic One-dimensional Burgers Equation

To investigate accuracy of the proposed DBFE method, its numerical results are compared with the Monte Carlo simulation. The Monte Carlo method for a stochastic one-dimensional Burgers equation is applied by solving the deterministic Burgers equation with initial condition defined using samples of . In this section, results are presented for a test case with . For each sample, explicit fourth order Runge-Kutta scheme is used for time evolution while the first order central KT scheme [20] (see equation (LABEL:num_flux)) is used in spatial dimension. The numerical solutions are obtained for and . Total 10000 samples are collected for the Monte Carlo simulation. Figure LABEL:ini_cond_mean(a) shows several realizations of initial conditions, where the mean is shown with emphasis, while respective solutions at time are shown in Figure LABEL:ini_cond_mean (b). Shock formation is observed for all the samples, however, shock location varies depending on the initial condition.

Figure 4.1: Solution of SPDE (LABEL:spde) using Monte Carlo method. Figure a) shows samples of initial condition; figure b) shows solution of samples corresponding to the initial conditions in a).

4.2 Application Results

Use the bi-orthogonal expansion of in (LABEL:spde_be) to obtain


The differential operator (LABEL:f_be) is used in (LABEL:num_flux_dbfe) to derive the DBFE governing equations (LABEL:mean_ev-LABEL:eyp_ev).

Initial condition for the mean is specified using (LABEL:ini_mean), while the initial condition for and are given by solution of the eigenvalue problem for the covariance function (LABEL:ini_cov). Note that the eigenvalues for the covariance function (LABEL:ini_cov) decreases rapidly with only first three eigenmodes dominant, thus, suffices for the accurate approximation of the uncertain initial condition. However, to account for a highly nonlinear nature of the Burgers equation, higher value of may be required. Figure LABEL:ini_eigen shows eigenvalues and first four eigenfunctions of the uncertain initial condition with and .

Figure 4.2: Eigenfield for uncertain initial condition. a) eigenvalues and b) first four eigenfunctions

4.2.1 Non-oscillatory Solution using the Proposed Numerical Scheme

Efficacy of the numerical scheme described in the section 3.1 to obtain non-oscillatory solutions for mean and eigenfunctions is investigated by implementing the DBFE method using the first order central KT scheme (LABEL:num_flux_dbfe), central difference scheme (terms involving , and are neglected in (LABEL:num_flux_dbfe)) and the first order central KT scheme in mean (terms involving and are retained only for the first term in (LABEL:num_flux_dbfe)).

Figure LABEL:comp_gibbs shows comparison of results for (a) mean and (b) first eigenfunction for . With central difference scheme, oscillations are observed in both mean and eigenfield. When central KT scheme is used only in mean field and central difference for eigenfield, oscillations in mean are comparatively reduced. However oscillations are still present in eigenfield, as can be observed for first eigenfunction. When full central KT scheme is used, oscillations in both mean and eigenfield are resolved. Note that only averaged statistics in terms of expected values and covariance functions affect time evolution of mean and eigenfields, hence, no special treatment is necessary for random field.

Figure 4.3: Figure shows comparison of numerical solution for a) mean and b) first eigenfunction, obtained by implementing the proposed DBFE method using the central difference scheme, first order KT scheme in mean and full first order central KT scheme. The results are obtained using the first three eigenmodes (N=3) in bi-orthogonal expansion.

4.2.2 Characteristics of DBFE Solutions

The characteristics of the numerical solution for DBFE before post-processing is investigated. In order to investigate effect of truncation of KL expansion, Burgers equation is solved for , and . To compare results with Monte Carlo simulations, is reconstructed for every sample using (LABEL:kl_expn_dbfe). Accuracy of the DBFE method is compared with Monte Carlo using relative error in norm, which is given by

Remark 1

Note that the -error in (LABEL:l1_error) is defined for a small domain , and thus is not a point error. In the present paper, the -error at a point is defined using the adjoining two cells.

Figure LABEL:comp_err shows the -error of the DBFE method relative to the Monte Carlo method. error of the order of is observed away from the shock location, however, near the shock location, error of the order of is observed that does not reduce significantly with increasing .

Figure 4.4: Figure shows error for DBFE method. The error is calculated against the Monte Carlo method using 10000 samples.

Figure LABEL:comp_eigen (a) shows the first eigenfunction at for . The eigenfunction shows periodic oscillations near shock, with number of modes increasing with . Figure LABEL:comp_eigen (b) shows evolution of variance of different modes over time for . Note that the modes with negligible variance at become dominant over time evolution. This effect can be significantly observed for the eigenmode, which has negligible variance at but becomes the second most dominant mode over the time evolution. Thus, it is necessary to use higher number of modes, even though the modes show negligible variance at .

Figure 4.5: Figure shows a) 1st eigenfunction for , and , and b) variance for all modes (N=7)

Figure LABEL:comp_conf shows comparison of the 90% confidence bound for the Monte Carlo and the DBEE method. Note that the confidence bound shows significant disagreement near shock location for all the test cases. The discrepancy in the confidence bound is due to the oscillations in the samples near the shock location, characterizing the Gibbs phenomenon for each sample.

Figure 4.6: Figure shows comparison of 90% confidence bound.

The key observations from the results presented here can be summarized as: a) numerical solution of the DBFE exhibits oscillations near the discontinuity, b) -error is maximum at the shock location, c) The amplitude of the maximum -error does not decrease by increasing the number of eigenfunctions, . These observations are characteristics of the Gibbs phenomenon, indicating the need for appropriate post-processing. Efficacy of the proposed post-processing to mitigate the effect of the Gibbs phenomenon is demonstrated in the following.

4.2.3 Mitigation of Gibbs Phenomenon using Post-processing

The Gegenbauer reprojection based approach, outlined in section 3.2, is used for post-processing to mitigate the effect of the Gibbs phenomenon. The bi-orthogonal expansion of is reprojected on order Gegenbauer polynomials, while total 7 coefficients are used in the expansion (LABEL:gegen_part_sum_pc). Integrals involving the Gegenbauer polynomials are numerically evaluated using Gauss-Gegenbauer quadrature with 100 nodes, whereas, the integrals involving the polynomial chaos basis are numerically solved using Monte Carlo integration. In the present paper, shock is assumed to be located at a point with highest slope where line is crossed. However, the proposed method can be implemented using any edge detection method for shock localization.

Resultant samples after post-processing using reprojection method are compared with samples without post-processing in Figure LABEL:post_process (the first six samples are shown in the figure). The respective Monte Carlo samples are also shown in the figure. Resolution of the Gibbs phenomenon can be observed for individual samples. For all the samples, oscillations are observed for approximations of obtained using DBFE method. Magnitude of oscillations is low for samples near mean and increases away from mean. The oscillations show three distinct behaviors: (1) when shock is located at , oscillations are observed on left side of shock location; (2) when shock is located at , oscillations are observed on right side of the shock location; and (3) when shock is located near , oscillations are observed on both the sides with low magnitude.

Figure 4.7: Figure shows comparison of DBFE method with Monte Carlo samples. Comparison without and with post-processing are shown for first six samples.

Figure LABEL:comp_gb a) shows comparison of 90% confidence bound of DBFE method after post-processing with Monte Carlo simulation. Using post-processing, agreement with Monte Carlo simulation is significantly improved. Figure LABEL:comp_gb b) shows error after post-processing. Comparing with Figure LABEL:comp_err, it can be observed that the error is reduced at all locations by the order of . At mean shock location, error is of the order of .

Figure 4.8: Figure shows a) confidence bounds using post-processing b) error

Figure LABEL:comp_moments shows comparison of first four moments of DBFE solution with the Monte Carlo simulation. The comparison is shown for both without and with post-processing cases. For the higher moments, agreement with the Monte Carlo method improves significantly after the post-processing. Mean of the solution obtained using the DBFE method without post-processing matches well with the Monte Carlo method (shown in Figure LABEL:comp_moments (a)). Note that the mean is not contaminated by the Gibbs phenomenon, resulting in the good match between DBFE method without post-processing and the Monte Carlo method, while, the subsequent post-processing retains the accuracy of the mean. Figure LABEL:comp_moments (b) shows the variance as a function of spatial location . The variance is low in the region away from the shock location, whereas, significant increase in the variance is observed at mean shock location. Good agreement between the Monte Carlo and the DBFE method is observed in the region away from the shock location, however, nontrivial difference is observed without post-processing near the mean shock location, with DBEF predicting higher variance than the Monte Carlo. The high variance in DBFE predictions results due to the oscillations near shock location characterizing the Gibbs phenomenon, which are mitigated using the post-processing, providing the close agreement with the Monte Carlo method. Figure LABEL:comp_moments also shows comparison for (c) skewness and (d) kurtosis. Skewness is zero in a region away from the shock location, indicating symmetric probability distribution. In a region near shock, skewness is negative in the left side of the mean shock location and positive on the right, with significantly high absolute value near the mean shock location . Similar results are obtained for kurtosis (Figure LABEL:comp_moments (d)) where high values are obtained near shock. It can be observed that DBFE method has predicted the trend correctly for skewness and kurtosis, while the close agreement with the Monte Carlo method is obtained after post-processing.

Figure 4.9: Figure shows comparison of first four moments for Monte Carlo and DBEF. a) mean b) variance c) skewness d) kurtosis

4.3 Computational Efficiency of the Proposed Method

The computational cost of the numerical implementation of the proposed DBFE method is compared with the Monte Carlo and the gPC method. See Pettersson et al. [35] for the governing equations of the gPC method applied to stochastic Burgers equation. The method is implemented with order gPC basis. Figure LABEL:conf_bound_gpc shows the 90% confidence bound for the gPC method, which is compared with the Monte Carlo method. Comparison for the mean is also shown in the figure. The mean obtained using the gPC method agrees well with the Monte Carlo method, however, the confidence bound shows oscillation near the mean shock location, demonstrating the existence of the Gibbs phenomenon for solution of the SPDE (LABEL:spde) using the gPC method. Note that the similar observations are reported earlier in the literature for solution of the stochastic Burgers equations using the gPC method [35].

Figure 4.10: Figure shows comparison of the confidence bound obtained using the gPC method with the Monte Carlo method. Comparison for mean is also shown in the figure.

To investigate the computational efficiency of the proposed method, the Monte Carlo method, the gPC method and the DBFE are implemented on a desktop computer with Intel Core i5 CPU. Total 10000 samples are used for the Monte Carlo method, while, order polynomial chaos basis is used for the DBFE and the gPC method. Time required for simulations is obtained using the FORTRAN intrinsic routine cpu_time. The comparison of CPU time is shown in Figure LABEL:comp_time. Both the gPC and the DBFE methods provide considerable computational speed-up over the Monte Carlo method, with the computational time requirement increasing with the number of eigenfunctions used. Even for , speed up of the order of 3 is obtained as compared to Monte Carlo method. The computational cost of the DBFE method is lower than the gPC method, with the ratio of the CPU time for gPC and the DBFE decreasing with increasing . Although the post-processing increases the computational cost of the DBFE method, the DBFE method still remains computationally efficient with the CPU time approaching that of the gPC with . Note that the computational cost of the proposed method depends on the complexity of the differential operator , however, the computational cost of the post-processing remains constant irrespective of the complexity of the SPDE.

Figure 4.11: Figure shows comparison of computational time for simulation. 10000 samples are used for the Monte Carlo method. order Hermite polynomials are used as basis for DBFE and gPC.

4.4 Choice of Gegenbauer Parameters

The proposed reprojection method requires specification of parameter of and maximum number of expansion terms . Total error in the Gegenbauer reprojection method eminates from error due to approximation of a function using a finite Gegenbauer series expansion and roundoff error in computation of Gegenbauer polynomials . Roundoff error in evaluation of increases with increasing and  [7], rendering very large and unsuitable.

To investigate the effect of choice of parameters of Gegenbauer polynomial, proposed post-processing method is applied using different values of and . Figure LABEL:tol_nk shows minimum value of for which error is greater than the tolerance for different . For given , limiting value of decreases with . In the present study, is used. note that since is linear on the both side of shock, is dominant for all the samples, while values of higher coefficients are very low. The higher coefficients control non-linearity near shock. For , highest order of Gegenbauer polynomial is quadratic in nature, which can not properly capture non-linearity near shock. Since the values of coefficients decrease rapidly, roundoff errors dominates for higher Gegenbauer polynomials resulting in resulting in oscillations near shock discontinuity. For the test case presented in this paper, choice of and provide trade-off between requirement for capturing non-linearity near shock and lower round-off error. Although the choice of parameters is heuristic in nature, a robust optimum parameter selection method can be developed based on already established techniques for deterministic solutions [7].

Figure 4.12: Figure shows limiting value of for given tolerance.

4.5 Robustness of the Method

Robustness of the proposed method is investigated for following covariance functions: 1) squared exponential , 2) triangular and 3) uniformly modulated . Figure LABEL:comp_err_sigma shows error for different values of when squared exponential covariance function is used to represent uncertainty in initial conditions. error increases exponentially with increasing , while, the error decreases after post-processing for all . For high , coefficients of higher modes of eigenfield become significant, necessitating the use of high . Thus, error due to truncation of spectral expansion at terms dominates the error for high , resulting in lesser improvement using post-processing.

Figure 4.13: Figure shows error as function of the variance .

Figure LABEL:conf_bound_tc shows comparison of 90% confidence bound obtained using MC and DBFE for squared exponential, triangular and uniformly modulated covariance function. Confidence bounds for DBFE method agrees well with MC, demonstrating robustness of the proposed method to resolve Gibbs phenomenon for different choices of covariance functions.

Figure 4.14: Figure shows comparison of confidence bound for different test cases.

5 Concluding Remarks

This paper have proposed a dynamic bi-orthogonality based spectral projection method for uncertainty quantification for systems with discontinuous solutions. The proposed approach is implemented in two steps: in the first step, input uncertainty is propagated to the system response using DBFE method, while, in the second step, effect of the Gibbs phenomenon is mitigated by reprojecting the mean and the eigenfield on the Gegenbauer polynomials. Efficacy of the proposed method have been investigated for solution of a one-dimensional stochastic Burgers equation. The numerical results presented in this paper have demonstrated the ability of the proposed post-processing method to mitigate the effect of the Gibbs phenomenon as discontinuities develop in the solution. Note that the proposed method does not require a-priory knowledge of the shock location, thus, a generic implementation for a variety of applications can be achieved by extending any numerical scheme to the DBFE, as demonstrated in this paper. The DBFE method is found to be computationally efficient than the gPC method, thus, the method becomes an attractive alternative to the gPC for solution of SPDEs.


  • [1] R.J. Adler and J.E. Taylor, Random Fields and Geometry, Springer, New York, 2007.
  • [2] R.H. Cameron and W.T. Martin, The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals, The Annals of Mathematics, 48 (1947), pp. 385–392.
  • [3] A. Chakrabarti and S.C. Martha, Approximate solutions of Fredholm integral equations of the second kind, Applied Mathematics and Computation, 211 (2009), pp. 459–466.
  • [4] T. Chantrasmi, A. Doostan, and G. Iaccarino, Pade-legendre approximants for uncertainty analysis with discontinuous response surfaces, Journal of Computational Physics, 228 (2009), pp. 7159–7180.