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 postprocessing 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 (DO) 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 postprocessing 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 and stochastic quasi-one-dimensional flow through a convergent-divergent nozzle.

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 [14]. Subsequently, the research community has invested significant efforts to develop various uncertainty quantification methodologies (see Ref. [5] 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.

Background.

For further discussion, consider a probability space (Ω,F,P) where ωΩ is a set of elementary events, F is associated σ-algebra, and P is a probability measure defined over F. A generic stochastic partial differential equation (SPDE) on (Ω,F,P) is given by 
u(x,t;ω)t=L[u(x,t;ω);ω]
(1)
where ω denotes stochasticity, u(x,t;ω) is a random variable, tT is time, xXRd is a spatial dimension, and L is an arbitrary differential operator of the form 
L[u(x,t;ω);ω]=f(u(x,t;ω);ω)x
(2)
and f(u(x,t;ω);ω) being the flux. The initial condition of the system is specified using a random field u(x,0;ω), while the boundary conditions are given by 
B(β,t;ω)=h(β,t;ω);   βX, ωΩ
(3)

where B is a differential operator.

Owing to the generality and ease of implementation, Monte Carlo methods are widely used for solution of Eq. (1). The method involves simulation of the system at large number of conditions sampled from input probability distribution. The efficiency of the Monte Carlo method can be improved by using variance reduction techniques like importance sampling [6] or sequential Monte Carlo methods [7], among others. Note that the Monte Carlo methods exhibit the accuracy proportional to Ns1/2 [6], where Ns is the number of samples collected. Thus, the Monte Carlo method typically requires a large number of samples for acceptable accuracy [6]. For the simulation of large-scale systems, the 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. Note that the dynamical response of the SPDE, u(x,t;ω), can be represented as a mapping 
u:T×X×FR
(4)
Under a generic condition of square integrability, u(x,t;ω) can be represented as a parametric field in a Hilbert space H(T×X×F) [8]. The SSP methods often use the tensor product representation of H(T×X×F) in appropriate Hilbert spaces H1 and H2. For H(T×X×F), possible choices of tensor product representation are [9] 
H(T×X×F)=H1(T×F)×H2(X)
(5)
 
=H1(T)×H2(X×F)
(6)
 
=H1(T×X)×H2(F)
(7)

Note that choice of the orthogonal basis is made based on the Hilbert spaces H1 and H2. A Karhunen–Loeve expansion [10] uses the tensor product (5), where u(x,t;ω) is spectrally represented in terms of orthogonal basis in H2(X) with temporally varying expansion coefficients in H1(T×F). A generalized polynomial chaos (gPC) expansion [11] is an example of tensor product (7) with orthogonal basis in H2(F) and the associated expansion coefficients in H1(T×X). This paper investigates a bi-orthogonal method, which uses the gPC basis in H2(F), while the dynamically orthogonal eigenfunction basis is used in H1(T×X).

The present state of the art for SSP methods is based on the Galerkin projection-based generalized polynomial chaos method introduced by Xiu and Karniadakis [11,12], which uses the representation of form (7). Using separability, H1(F) is decomposed in orthogonal polynomial chaos basis ψk(ω), where ψk(ω) belongs to the Askey scheme of orthogonal polynomials [13] with weight function given by probability distribution function of an uncertain input. Thus, using gPC, the dynamical response u(x,t;ω) is spectrally represented as 
u(x,t;ω)=i=1Pui(x,t)ψi(ω)
(8)

where ui(x,t) are dynamical expansion coefficients.2 See Mathelin et al. [14] and Najm [15] for extensive review of the research work in this field.

Motivation.

The present work is motivated by the need to quantify uncertainty in simulation of systems with discontinuous solutions. An example of such a discontinuity is shocks that evolve during transonic simulations. In particular, the paper deals with stochastic systems defined using the hyperbolic conservation laws that invariably result in a discontinuous solution irrespective of the initial conditions. The simulation of such system has attracted significant interest by the research community and the literature is rich with methods for resolution of discontinuities in deterministic simulations (see Ref. [16] and references therein for a review of the methods). However, the problem of the development of numerical methods for solution of stochastic hyperbolic partial differential equations is not yet resolved and is of current interest to the research community [17].

Considering the cost of the Monte Carlo methods, the present paper focuses on SSP methods for the solution of SPDE in Eq. (1). The use of spectral methods for hyperbolic SPDEs is known to pose the following significant numerical challenges [18]

  1. (1)

    Numerical solutions using finite difference schemes lead to spurious oscillations in the expansion coefficients and/or orthogonal basis.

  2. (2)

    The spectral expansion of the solution of hyperbolic SPDEs lead to Gibbs phenomenon as discontinuities develop in solution [18]. The Gibbs phenomenon for application of SSP methods to hyperbolic SPDEs with discontinuous solution is characterized by

    • (a)

      Slow convergence of the spectral approximation of the solution at points away from the discontinuity. The convergence is O(1/N), where N is the number of expansion terms retained in the spectral approximation.

    • (b)

      O(1) oscillations are observed in the solution near discontinuities that do not decrease with increasing N.

The Gibbs phenomenon is discussed in the deterministic settings for Fourier analysis by Hewitt and Hewitt [19] and for more generic orthogonal polynomials by Gottlieb and Shu [20].

Proposed Method.

The existence and resolution of Gibbs phenomenon in spectral expansion of a deterministic function using a finite number of orthogonal bases is extensively investigated in the literature [1921]. The spectral expansion of a function using global orthogonal basis is contaminated by the presence of discontinuities. However, the expansion coefficients contain enough information to recover the function with high accuracy using appropriate postprocessing. Gibbs phenomenon can be completely resolved by reprojecting the partial sum on Gibbs complementary basis [20]. Gottlieb et al. [2023] have shown that exponentially convergent approximation can be obtained at point values in a subinterval of analyticity by reprojecting the partial sum on orthogonal Gegenbauer polynomial basis. The postprocessing method can be extended to stochastic functions, provided, the spectral expansion is available in terms of orthogonal basis in spatial dimensions.

Note that the generalized polynomial chaos (gPC) expansion (8) uses an orthogonal basis in H2(F) while corresponding coefficients ûi(x,t) are functions in H1(T×X); however, the discontinuity resides in the spatial dimension. Thus, enough information is not available in ûi(x,t) to resolve Gibbs phenomenon using the Gegenbauer reprojection method. The resolution of the Gibbs phenomenon proposed in this paper is based on the observation that exponential accuracy can be recovered if u(x,t;ω) is expanded in terms of orthogonal bases in the spatial dimension, and subsequently reprojected on an approximate Gibbs complementary basis. Consider a generic Karhunen–Loeve expansion of u(x,t;ω) 
u(x,t;ω)=u¯(x,t)+i=1NYi(t;ω)ui(x,t)
(9)

where u¯(x,t) is the mean, ui(x,t) are eigenfunctions which form a complete orthogonal basis in H1(X) at a given time step t, while Yi(t;ω) are independent zero-mean random variables. Sapsis and Lermusiaux [24] have proposed the dynamic orthogonality (DO) condition, which constrains (ui(x,t))/(t) to be orthogonal to uj(x,t), to ensure orthogonality of eigenfunctions at all time steps.

Tagade and Choi [25] 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 with discontinuous solution using DBFE, which exploits orthogonal decomposition of the spatial dimension to develop a postprocessing approach for mitigation of the Gibbs phenomenon. The DBFE method represents solution of a stochastic system at each time step as the decomposition in mean and random field. The random field is represented as convolution of Hilbert spaces in stochastic and spatial dimensions. Using polynomial chaos expansion of Yi(t;ω) in Eq. (9), bi-orthogonal expansion of u(x,t;ω) is given by 
u(x,t;ω)=u¯(x,t)+i=1Np=1PYpi(t)ψp(ω)ui(x,t)
(10)

where ψp(ω) are orthogonal polynomials. The DO condition is used to derive closed-form solution of evolution equations for u¯(x,t),ui(x,t), and Ypi(t). However, the resultant approximation of u(x,t;ω) in Eq. (10) using the solution of u¯(x,t),ui(x,t) and Ypi(t) exhibits the Gibbs phenomenon. A Gegenbauer reprojection-based postprocessing method is proposed to mitigate the Gibbs phenomenon. During the postprocessing step, eigenfunctions ui(x,t) are reprojected on the Gegenbauer basis. Note that the Bi-orthogonal expansion is explored by Cheng et al. [26], where the independence of expansion coefficients is used to obtain Bi-orthogonality. Choi et al. [27] combine the gPC expansion with the DO method, where the gPC method is used during the initial time evolution, while the DO expansion is used subsequently when the uncertainty in the solution is high enough. Key contribution of the present paper is a Gegenbauer reprojection method for mitigation of the Gibbs Phenomenon.

There have been recent research efforts that focused on addressing the issue of Gibbs phenomenon in gPC settings. Wan and Karniadakis [28] have proposed a multi-element gPC (ME-gPC) method to mitigate the effect of Gibbs phenomenon (also see Lin et al. [29]). 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. The ME-gPC method uses localized expansion in the subdomains to mitigate the Gibbs phenomenon. However, the computational cost of the ME-gPC method rises quickly as the number of subdomains increases. Poette et al. [30] have proposed an entropy-based intrusive polynomial moment method for uncertainty propagation in nonlinear systems with discontinuous solutions. This method exhibited better computationally efficient than the ME-gPC, but the requirement for minimization of the entropy at each time step results in the computational cost higher than the gPC method. 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 more efficient than the ME-gPC; however, the requirement for minimization of the entropy at each time step results in a computational cost higher than the gPC method. Tryoen et al. [31] use a Roe-type solver for intrusive Galerkin projection of uncertain hyperbolic systems. The solver uses the Dubois and Mehlman entropy corrector that avoids certain entropy violating solutions. Sargsyan et al. [17] 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 analyticity, while intrusive Galerkin projection is used to propagate the uncertainty. In an alternate non gPC setting, Chantrasmi et al. [32] have proposed a Pade–Legendre interpolant-based approach, where simulation runs are obtained at predefined Gauss–Legendre quadrature nodes and resultant discontinuous system response is reconstructed using the Pade interpolation. Despite this recent progress, the resolution of Gibbs phenomenon for stochastic simulation of systems with discontinuous solution using SSP methods has not been fully resolved and is an area of current research interest [17]. In addition, in the methods based on the gPC setting, the computational cost grows quickly as the stochastic dimension increase.

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

Dynamically Bi-Orthogonal Field Equations (DBFE)

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

Definition 1. Inner product onH1(T×X)for a given t is defined as 
u(x,t;ω),v(x,t;ω)X=Xu(x,t;ω)v(x,t;ω)dx
(11)
Similarly, inner product onH2(F)is defined as 
u(x,t;ω),v(x,t;ω)Ω=Ωu(x,t;ω)v(x,t;ω)dP(ω)
(12)
The expectation operator onH2(F)is defined as 
Eω[u(x,t;ω)]=Ωu(x,t;ω)dP(ω)
(13)
Further using definition of the expectation operator, the covariance is defined as 
Cu,v=Eω[u(x,t;ω)v(x,t;ω)]
(14)
The DBFE method uses the bi-orthogonal expansion (10) of u(x,t;ω) in Eq. (1), which results in 
u¯(x,t)t+i=1Np=1Pψp(ω)tYpi(t)ui(x,t)=L[u(x,t;ω);ω]
(15)
Note that Eq. (15) is a set of N×P+1 coupled partial differential equations (PDEs), thus, independent evolution equations for the mean u¯(x,t), the eigenfunctions ui(x,t) and the corresponding coefficients Ypi(t) cannot be obtained concurrently using the SPDE in Eq. (1). Well-posed evolution equations for the unknown quantities can be derived by imposing the additional dynamic orthogonality (DO) condition [24]. The DO condition is specified as [24] 
ui(x,t),uj(x,t)tX=0,  i, j=1,...,N
(16)

where N is the number of eigenfunctions retained in the expansion (10). 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 [24].

Using the DO condition and bi-orthogonal expansion of u(x,t;ω), SPDE (1) can be reformulated into a set of N + 1 PDEs and N × P ordinary differential equations (ODEs) as follows.

Proposition 1. Using the DO condition, dynamic evolution equations ofu¯(x,t),ui(x,t)andYpi(t)are given by 
u¯(x,t)t=Eω[L[u(x,t,ω);ω]]
(17)
 
ui(x,t)t=j=1NCYiYj1{Eω[L[u(x,t,ω);ω]Yj(t;ω)]k=1NEω[L[u(x,t,ω);ω]Yj(t;ω)],uk(x,t)X}
(18)
 
dYpi(t)dt=1ψp2(ω)ΩL[u(x,t;ω);ω]Eω[L[u(x,t;ω);ω]],ui(x,t)X,ψp(ω)Ω
(19)

Proof. Proof of Eqs. (17) and (18) is given by Sapsis and Lermusiaux [24], which is briefly presented here for completeness, while Eq. (19) is derived here by introducing the bi-orthogonality.

Proof of Eqs.(17)and(18). Use a generic Karhunen–Loeve (KL) expansion (9) in Eq. (1) to obtain 
u¯(x,t)t+i=1Nui(x,t)dYi(t;ω)dt+i=1NYi(t;ω)ui(x,t)t=L[u(x,t;ω);ω]
(20)

Application of the expectation operator to Eq. (20) gives (17).

Multiply (20) by Yj(t;ω) and apply the expectation operator to obtain 
i=1NCdYi(t)dtYj(t)ui(x,t)+i=1NCYi(t)Yj(t)ui(x,t)t=Eω[L[u(x,t;ω);ω]Yj(t;ω)]
(21)
Multiply (21) by uk(x,t), take inner product, and apply the DO condition to get 
CdYk(t)dtYj(t)=Eω[L[u(x,t;ω);ω]Yj(t;ω)],uk(x,t)X
(22)
Use Eq. (22) in Eq. (21) to obtain 
i=1NCYi(t)Yj(t)ui(x,t)t=Eω[L[u(x,t,ω);ω]Yj(t;ω)]k=1NEω[L[u(x,t,ω);ω]Yj(t;ω)],uk(x,t)Xuk(x,t)
(23)
which can be written in the matrix form as 
U=Σ1D
(24)

where Σ is the covariance matrix with ijth element Σij=CYiYj.

Multiply both sides of Eq. (20) by uj(x,t), take inner product, and use the DO condition and orthonormality of eigenfunctions to get 
u¯(x,t)t,uj(x,t)X+dYj(t;ω)dt=L[u(x,t;ω);ω],uj(x,t)X
(25)
which on application of expectation operator gives 
u¯(x,t)t,uj(x,t)X=Eω[L[u(x,t;ω);ω],uj(x,t)X]
(26)
Using Eq. (26) in Eq. (25) gives 
dYi(t;ω)dt=[L[(x,t;ω);ω]Eω[L[u(x,t;ω);ω]]],ui(x,t)X
(27)
Proof of Eq.(19). Basic conditions of KL expansion ensure that Yi(t;ω) are square integrable random variables; thus, according to Cameron and Martin theorem [33], Yi(t;ω) can be approximated to arbitrary accuracy using spectral expansion in terms of orthogonal basis in H2(F). Hence, Yi(t;ω) can be spectrally represented using P terms as [11] 
Yi(t;ω)=p=1PYpi(t)ψp(ω)
(28)
where ψp(ω) are orthogonal basis in H2(F). Differentiating (28) with respect to t gives 
dYi(t;ω)dt=p=1PdYpi(t)dtψp(ω)
(29)
which on Galerkin projection provide 
dYpi(t)dt=dYi(t;ω)dt,ψp(ω)Ωψp2(ω)Ω
(30)

where ψp2(ω)Ω denotes double product. Having (29) and (30) in Eq. (27) results in Eq. (19).

Initial and Boundary Conditions for DBFE.

The initial conditions for DBFE are defined by specifying the mean u¯(x,0), the eigenfunctions, ui(x,0) and the coefficients Ypi(0), which are given by considering the bi-orthogonal expansion of the initial condition of SPDE (1) as 
u(x,0;ω)=u¯(x,0)+i=1Np=1PYpi(0)ui(x,0)ψp(ω)
(31)
Applying the expectation operator, initial condition for the mean field is given by 
u¯(x,0)=Eω[u(x,0;ω)]
(32)
The initial condition for eigenfunctions, ui(x,0), is given by solution of the Fredholms' integral equation of the second kind [34] 
XCu(x1,0),u(x2,0)ui(x1,0)dx1=λi2ui(x2,0)
(33)
where Cu(x1,0),u(x2,0) is the covariance function of u(x,0;ω) and λi are the eigenvalues. A Galerkin projection-based method is used to solve (33) numerically [35]. For a Gaussian process, all the coefficients Ypi(0) are zero except 
Y2i(0)=λi
(34)
thus defining Yi(0;ω) as normal variables with variance λi. For a generic non-Gaussian case, the initial expansion coefficients Ypi(0) are given by 
Ypi(0)=u(x,0;ω)u¯(x,0),ui(x,0)X,ψp(ω)Ωψp2(ω)Ω
(35)
Boundary conditions for DBFE are given by using the generic KL expansion of h(β,t;ω), which specifies the boundary conditions for the SPDE in Eq. (1), as 
h(β,t;ω)=h¯(β,t)+i=1NYi(t;ω)ui(β,t)
(36)
Applying the expectation operator on Eq. (36), boundary condition for the mean field is given by 
B[u¯(β,t;ω)]βX=h¯(β,t)
(37)
Multiply (37) by Yj(t;ω) and apply the expectation operator 
Eω[h(β,t;ω)Yj(t;ω)]=i=1NCYi(t)Yj(t)ui(β,t)
(38)
which can be specified in a matrix form as 
u=Σ1E
(39)

where Σ is a covariance matrix.

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 u¯(x,t) and the eigenfunctions (using Eqs. (17) and (18)) results in the development of spurious oscillations as discontinuities evolve in the solution; (b) the resultant bi-orthogonal expansion (10) 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 nonoscillatory 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.

Numerical Scheme for Nonoscillatory Solution.

Due to discontinuities, numerical solution of Eqs. (17)(19) is expected to contain spurious oscillations with reduced accuracy in u¯(x,t) and ui(x,t). Total variation diminishing (TVD) finite difference schemes are widely used in the literature to obtain nonoscillatory accurate weak solutions to hyperbolic partial differential equations [16]. This paper demonstrates extension of a TVD scheme for numerical solution of Eqs. (17)(19); however, note that other numerical schemes, like essentially nonoscillatory schemes, Reimann solvers, Godunov approach, etc. (see Ref. [16] for details), can similarly be extended for DBFE.

For notational convenience, this subsection defines ujt(ω) as an approximate value of u(x=xj,t;ω) at the grid point xj=jΔx and present time t. Consider a numerical scheme applied to SPDE (1) for a given ω as 
ujt(ω)t=1Δx(f̂j+12(ω)f̂j12(ω))
(40)
where f̂j+12(ω) is an interpolated flux inside a cell [xj,xj+1]. A (2k+1) point numerical scheme is defined by using k points on either side of xj for interpolation, thus 
f̂j+12(ω)=f(ujk+1t(ω),...,uj+kt(ω))
(41)
such that [16] 
f̂(u,...,u)=f(u)
(42)
f̂j12(ω) is also defined in a similar manner. The TVD scheme (40) can be extended to DBFE by using 
L[u(x,t;ω);ω]=1Δx(f̂j+12(ω)f̂j12(ω))
(43)

for numerical solution of Eqs. (17)(19). The resultant numerical scheme retains nonoscillatory property in u¯(x,t) and ui(x,t).

In the present paper, the DBFE method is implemented using the first-order central scheme proposed by Kurganov and Tadmor [36] (central KT scheme), for which numerical flux is defined as 
f̂j+12(ω)=12[f(uj+1t(ω))+f(ujt(ω))]+12[aj+12(ω)(uj+1t(ω)ujt(ω))]
(44)
where aj+12(ω) is the maximum local speed at xj. The gPC expansion of aj+12(ω) is given by 
aj+12(ω)=p=1Pâpj+12ψp(ω)
(45)
where âpj+12 are gPC expansion coefficients. Use Eq. (45) in numerical flux (44) to obtain 
Eω[f̂j+12(ω)]=12{Eω[f(uj+1t(ω))]+Eω[f(ujt(ω))]}+12{â1j+12[u¯(xj+1,t)u¯(xj,t)]+i=1NM1i,j+12[ui(xj+1,t)ui(xj,t)]}Eω[f̂j+12(ω)Yk(t;ω)]=12{Eω[f(uj+1t(ω))Yk(t;ω)]+Eω[f(ujt(ω))Yk(t;ω)]}+12{(M1k,j+12[u¯(xj+1,t)u¯(xj,t)]+i=1NT1i,k,j+12[ui(xj+1,t)ui(xj,t)]}f̂j+12(ω)Eω[f̂j+12(ω)]=12{f(uj+1t(ω))Eω[f(uj+1t(ω))]+f(ujt(ω))Eω[f(ujt(ω))]}+12{p=2Pâpj+12[u¯(xj+1,t)u¯(xj,t)]+i=1NT1i,k,j+12[ui(xj+1,t)ui(xj,t)]}
(46)

where Mi,j+12=Yi(t;ω)aj+12(ω) and Ti,k,j+12=Yi(t;ω)Yk(t;ω)aj+12(ω).

The numerical flux (46) is used in Eq. (40) to obtain the first-order central difference scheme for numerical solution of Eqs. (17)(19). Fourth-order Runge–Kutta method is used for time integration in Eqs. (17) and (18), while Euler method is used to solve ODE (19). Gauss–Legendre quadrature is used to calculate inner products in spatial dimensions, while Monte Carlo integration is used for inner products in stochastic dimensions.

Postprocessing Using Gegenbauer Polynomial.

Let the solution of the SPDE in Eq. (1) is obtained using the DBFE method till time t = T such that the discontinuities are developed. For a given ω, let [a(ω),b(ω)] denote an interval of analyticity for u(x,T;ω), where a random variable ζ(x;ω) can be defined such that 1ζ(x;ω)1. Gibbs phenomenon can be resolved in the interval [a(ω),b(ω)] by reprojecting the partial sum (10) on a suitable Gibbs complementary polynomial basis [21]. In the present paper, the partial sum is reprojected on a Gegenbauer polynomial, which is defined for a ζ(x;ω) as follows.

Definition 2. The Gegenbauer polynomialCnλ(ζ(x;ω))is a polynomial of degree n which is orthogonal over[1,1]with weight function(1ζ2(x;ω))λ12forλ0. The orthogonality is given by [21] 
11(1ζ2(x;ω))λ12Cnλ(ζ(x;ω))Cmλ(ζ(x;ω))dζ=hnλδm,n
(47)
where 
hnλ=πCnλ(1)Γ(λ+12)Γ(λ)(m+λ)
(48)
with 
Cnλ(1)=Γ(n+2λ)m!Γ(2λ)
(49)
Cnλ(ζ(x;ω))can be estimated using a recurrence relationship 
(n+1)Cn+1λ(ζ(x;ω))=2(λ+n)ζ(x;ω)Cnλ(ζ(x;ω))(2λ+n1)Cn1λ(ζ(x;ω))
(50)
Using orthogonality of Gegenbauer polynomial, u(x,T;ω) can be represented in terms of exponentially convergent Gegenbauer expansion series truncated at M terms as 
u(x,T;ω)=l=1Mûl(ω)Clλ(ζ(x;ω))
(51)
where ζ(x;ω)=xδ(ω)/ε(ω); ε(ω)=b(ω)a(ω)/2 and δ(ω)=b(ω)+a(ω)/2. Coefficients ûl(ω) are given by 
ûl(ω)=1hlλ11(1ζ2(x;ω))λ12Clλ(ζ(x;ω))u(x(ζ),T;ω)dζ
(52)
However, u(x(ζ),T;ω) is not available for calculation of ûl(ω). Using bi-orthogonal expansion (10) in Eq. (52), approximate Gegenbauer expansion coefficients ĝl(ω) are given by 
ĝl(ω)=1hlλ[11(1ζ2(x;ω))λ12Clλ(ζ(x;ω))u¯(x,T)dζ+i=1Np=1P(Ypi(T)ψp(ω)×11(1ζ2(x;ω))λ12Clλ(ζ(x;ω))ui(x,T)dζ)]
(53)
Note that (53) uses a bi-orthogonal expansion of u(x,T;ω); thus, ĝl(ω) is an approximation of ûl(ω). ĝl(ω) can be expanded in gPC basis as 
ĝl(ω)=p=1Pĝplψp(ω)
(54)
where 
ĝpl=ĝl(ω),ψp(ω)Ωψp2(ω)Ω
(55)
Coefficients ĝpl are calculated using Gaussian quadrature by evaluating (53) at collocation points, while the integrals in Eq. (53) are calculated using Gauss–Gegenbauer quadrature. Using ĝpl,u(x,T;ω) is approximated as 
u(x,T;ω)=l=1Mp=1Pĝplψp(ω)Clλ(ζ(x;ω))
(56)

The partial sum (56) approximates u(x,T;ω) with exponential accuracy in the subinterval [a(ω),b(ω)]. Note that determination of domain [a(ω),b(ω)] requires location of discontinuity. An edge detection method [37,38] can be used to locate points of discontinuities.

Note that the proposed method is implemented in two steps. The first step involves time evolution of the mean, eigenvalues, and eigenfunctions of the DBFE expansion. The proposed extension of the TVD scheme ensures unphysical conditions are not developed during this time evolution. However, the presence of the discontinuity in the solution results in poor quantification of uncertainty due to Gibbs phenomenon. Reprojection on the Gegenbauer polynomials recovers accuracy of this uncertainty quantification. The numerical implementation of this proposed approach is summarized in the Algorithm 1.

Algorithm 1 Numerical Implementation

1: procedure Time Evolution of DBFE Expansion

2:  Obtain DBFE expansion of initial and Boundary condition.

3:  Use numerical flux given by Eq. (46) with the Runge–Kutta method to obtain time evolution of mean, eigenfunctions and gPC expansion coefficients.

4:   At an arbitrary time instance T, use the DBFE expansion to obtain necessary statistics of the solution.

5: procedure Reprojection on Gegenbauer Polynomials

6:   At an arbitrary time instance T, use expansion (56) to reproject the solution on Gegenbauer polynomials.

7:     Obtain the DBFE expansion from the reprojected solution and restart the time evolution.

Numerical Examples

Stochastic 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 [39], while there have been recent investigation on solution of stochastic Burgers equation [30,40]. Coupled with the ease of implementation, stochastic Burgers equation provides 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 
t(u(x,t;ω))=x(u2(x,t;ω)2)
(57)
where x[1,1]. In the present paper, the stochastic Burgers equation is investigated for uncertain initial condition specified using a Gaussian process with mean 
u¯(x,0)=ubtan1(xx0)
(58)
where ub and x0 are user-defined constants, and the covariance function 
Cu(x1,0),u(x2,0)=σ2exp(λ|x1x2|)
(59)
where σ2 is variance and λ is the correlation length. The boundary condition is specified using 
u(x,0;ω)=u(β,0;ω)
(60)

Figure 1 shows deterministic solution of the Burgers equation at different time steps for mean initial condition with ub = 0 and x0=0. For the deterministic numerical solution, explicit fourth-order Runge–Kutta scheme is used for time evolution while the first-order central Kurganov and Tadmor scheme [36] (see Eq. (44)) is used in spatial dimension. One-dimensional grid is used with Δx=0.01 while time step of Δt=104 is used. With time, left state moves toward increasing x, while at time t*, when the characteristic lines first cross, solution has infinite slope at x = 0, while the solution becomes discontinuous at x = 0 for t>t*. Final solution is shown at time t=1.0s.

Solution of Stochastic Burgers Equation Using DBFE.

Use the bi-orthogonal expansion of u(x,t;ω) 
u(x,t;ω)=u¯(x,t)+i=1Np=1PYpi(t)ψp(ω)ui(x,t)
(61)
in Eq. (57) to obtain 
L[u(x,t;ω);ω]=xu¯2(x,t)2i=1Np=1PYpi(t)ψp(ω)x(u¯(x,t)ui(x,t))i=1Nj=1Np=1Pq=1PYpi(t)Yqj(t)ψp(ω)ψq(ω)x(ui(x,t)uj(x,t)2)

The differential operator (62) is used in Eq. (46) to derive the DBFE governing equations (17)(19). The field equations are spatially discretized with Δx=0.01, while the fourth-order explicit Runge–Kutta method with time step Δt=104 is used for the time integration.

Initial condition for the mean u¯(x,t) is specified using Eq. (58), while the initial condition for ui(x,t) is given by eigenfunctions of the covariance function (59). The initial conditions for expansion coefficients Ypi(t) are specified as per Eq. (34). Note that the eigenvalues for the covariance function (59) decrease rapidly with only first three eigenmodes dominant [10]; thus, N = 3 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 N may be required. As can be inferred from Fig. 1, solution of stochastic Burgers equation is expected to evolve discontinuities. Effect of the Gibbs phenomenon for the resultant DBFE solution is mitigated by (a) obtaining nonoscillatory solution for mean and eigenfield and (b) using postprocessing on the resultant spectral expansion.

Figure 2 investigates efficacy of the proposed central KT scheme (see Sec. 3.1) to obtain nonoscillatory solution using the DBFE method for the eigenfield. The figure compares solution of the eigenfield obtained using three different numerical implementations of the DBFE method: (a) central difference scheme, where terms involving â, M, and T are neglected in Eq. (46) (indicated as No Treatment in Fig. 2), (b) central KT scheme in mean that neglects terms involving â and M in Eq. (46) (Mean Treatment in Fig. 2), and (c) central KT scheme (All Treatment in Fig. 2). With central difference scheme, oscillations are observed in the eigenfield. When central KT scheme is used only in mean field and central difference for eigenfield, oscillations are still present in eigenfield, though the amplitude is reduced, as can be observed for the first eigenfunction. When full central KT scheme is used, oscillations in eigenfunction are resolved.

To investigate the accuracy of the proposed method, the DBFE solution of the stochastic Burgers equation is compared with the Monte Carlo method using 10,000 samples. Figure 3(a) shows comparison of the 90% confidence bound for Monte Carlo and the DBFE method. 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. Amplitude of the oscillations is negligible for the samples near mean, which increases with the distance of the samples from the mean. Note that the distance of a sample from the mean is characterized by expansion coefficients Ypi(t) and the associated gPC basis ψp(ω). For samples close to the mean, reconstruction of u(x,t;ω) through the bi-orthogonal expansion (10) is dominated by u¯(x,t), for which nonoscillatory solution is obtained using the numerical scheme (46). However, for samples away from the mean, effect of the Gibbs phenomenon becomes dominant through eigenfunctions ui(x,t). Figure 3(b) shows the L1-error of the DBFE method relative to the Monte Carlo method. L1(ω) error of the order of 102 is observed away from the shock location; however, near the shock location, error of the order of 101 is observed. Also, it is observed that the error does not reduce significantly with increasing N. From Figs. 3(a) and 3(b), it may be concluded that the L1-error is maximum at the shock location; and the amplitude of the maximum L1-error does not decrease by increasing the number of eigenfunctions, N. These observations are characteristics of the Gibbs phenomenon, indicating the need for appropriate postprocessing.

The Gegenbauer reprojection-based approach, outlined in Sec. 3.2, is used for postprocessing to mitigate the effect of the Gibbs phenomenon. The bi-orthogonal expansion of u(x,t;ω) is reprojected on seventh-order Gegenbauer polynomials, while total seven coefficients are used in the expansion (56). 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 u = 0 line is crossed. However, the proposed method can be implemented using any edge detection method for shock localization.

Figure 3(c) shows comparison of 90% confidence bound of DBFE method after postprocessing with Monte Carlo simulation. Using postprocessing, agreement with Monte Carlo simulation is significantly improved. Figure 3(d) shows L1 error after postprocessing. Comparing with Fig. 3(b), it can be observed that the error is reduced at all locations by the order of 102. At mean shock location, L1 error is of the order of 103.

Figure 4 shows comparison of the mean and variance of the DBFE solution with the Monte Carlo simulation. The comparison is shown for both without and with postprocessing cases. For the higher moments, agreement with the Monte Carlo method improves significantly after the postprocessing. Mean of the solution obtained using the DBFE method without postprocessing matches well with the Monte Carlo method (shown in Fig. 4(a)). Note that the mean is not contaminated by the Gibbs phenomenon, resulting in the good match between DBFE method without postprocessing and the Monte Carlo method, while the subsequent postprocessing retains the accuracy of the mean. Figure 4(b) shows the variance as a function of spatial location x. 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 postprocessing near the mean shock location, with DBFE 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 postprocessing, providing the close agreement with the Monte Carlo method.

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 Ref. [40] for details of the solution of stochastic Burgers equation using the gPC method). 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 10,000 samples are used for the Monte Carlo method, while third-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 Fig. 5. 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 N = 7, 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 N. Although the postprocessing 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 N. Note that the computational cost of the proposed method depends on the complexity of the differential operator L[u(x,t;ω);ω]; however, the computational cost of the postprocessing remains constant irrespective of the complexity of the SPDE. The proposed DBFE method with the postprocessing provides the solution to SPDE with the accuracy comparable to the Monte Carlo method at a significantly lower computational cost.

Quasi-One-Dimensional Nozzle Flow.

Solution of stochastic quasi-one-dimensional nozzle flow in the supersonic regime using the gPC method is investigated in the literature [41]; however, the gPC method fails for the subsonic-supersonic regime due to appearance of Gibbs phenomenon. This section investigates the proposed DBFE method for simulation of a quasi-one-dimensional nozzle flow with uncertain nozzle area in the subsonic-supersonic regime. The governing equations for the quasi-one-dimensional nozzle flow are given by the Euler system for compressible flow as 
ρAt+ρvAx=0ρvAt+ρv2Ax=ApxρEAt+(ρE+p)uAx=0
(62)

where A is nozzle area, ρ is density, v is velocity, E is total energy, and p is the pressure. For notational convenience, define u1=ρA,u2=ρvA, and u3=ρEA in the following.

Initial and Boundary Conditions.

In the present section, the proposed method is investigated for the uncertain nozzle area, while the initial conditions for ρ, v, and p are deterministic. The uncertainty in the nozzle area is specified using 
A(x;ω)G(A¯(x),ΣA)
(63)
where G(A¯(x),ΣA) represents a Gaussian process with mean A¯(x) and covariance function ΣA. In the present paper, squared exponential covariance function is used, which is given by 
ΣA(x1,x2)=σA2exp(x1x2)2
(64)
where σA2 is variance. To initialize the DBFE solution, use the Karhunen–Loeve expansion of A(x;ω) to obtain 
A(x;ω)=A¯(x)+i=1NAiai(x)
(65)
where ai(x) are eigenfunctions. Let ρini(x) be the initial condition for the density. Use the Karhunen–Loeve expansion 65 with ρini(x) to obtain 
u1(x)=ρini(x)A¯(x)+i=1NAiρini(x)ai(x)
(66)
Thus, initial condition for the mean of u1(x) is given by 
u1¯(x)=ρini(x)A¯(x)
(67)
Note that the covariance function for u1(x) is 
Σu1(x1,x2)=i=1Nj=1Nρini(x1)ρini(x2)ai(x1)aj(x2)ΣA(x1,x2)
(68)
Thus using Eq. (68), initial condition for the eigenfunction basis of u1(x),u1,i(x) can be obtained. The initial condition for the gPC basis of the DBFE expansion, then, can be obtained using 
Yp1,i(0)=u1(x,0;ω)u¯1(x,0;ω),u1,i(x,0),ψp(ω)
(69)

The DBFE expansion for u2(x,0;ω) and u3(x,0;ω) is also initialized similarly.

Deterministic boundary conditions are specified at both the nozzle inlet and exit conditions. At the inlet boundary, density and temperature are specified, while the velocity is linearly interpolated from the adjacent grid points. Static pressure is specified at the exit boundary condition, whereas density and velocity are extrapolated from the adjacent grid points. Note that due to uncertain nozzle area, u(B,t;ω) is uncertain despite the deterministic boundary conditions.

Numerical Results.

All the numerical results presented in this section are obtained using a first-order KT scheme in spatial dimension and fourth-order Runge–Kutta method for time integration, while spatiotemporal discretization is used with Δx=0.01 and Δt=105. The results are shown after total time of T=3s. Figure 6(a) shows the mean nozzle area profile, while Fig. 6(b) shows corresponding deterministic solution. Existence of the shock discontinuity in the divergent region of the nozzle can be observed, resulting in the subsonic flow at the nozzle exit.

Uncertainty in the nozzle area is specified using a Gaussian process defined in Eq. (63) with the variance σA2=0.1. The uncertain nozzle area is represented using a Karhunen–Loeve expansion with the mean nozzle profile as shown in Fig. 6(a), while N = 3 eigenfunctions are used for the spectral expansion. Resultant initial conditions for the DBFE solution are obtained as per the discussion in Sec. 4.3. Figure 7 shows the eigenfunctions at time t = 0 (i.e., initial condition) and t=3s (i.e., final solution). First eigenfunction is changed significantly for all u1(x,t;ω),u2(x,t;ω), and u3(x,t;ω); however, second and third eigenfunctions are close to the initial conditions for u1(x,t;ω) and u3(x,t;ω). However, note that all the eigenfunctions have evolved significantly over time for u2(x,t;ω). From the governing equations for quasi-one-dimensional nozzle flow, it can be noted that the uncertain nozzle area profile primarily affects the momentum equation; thus, significant evolution in all the eigenfunctions is observed for u2(x,t;ω).

To investigate efficacy of the proposed method, solution of stochastic quasi-one-dimensional nozzle flow is obtained using the Monte Carlo method. Relevant statistics is estimated by collecting 10,000 samples. Figure 8 shows comparison of mean velocity and mean pressure obtained using the Monte Carlo and the proposed method. As can be observed from Fig. 8(a), mean velocity and pressure obtained using the DBFE match closely with the Monte Carlo method. The relative L1-error in mean is shown in Fig. 8(b). The relative L1-error is of the order of 103, with the maximum error near the mean shock location for both the velocity and pressure. Note that for the reason of brevity, results only for velocity and pressure are shown in this section, though other parameters show similar trend.

Figure 9 shows comparison of variance obtained using the Monte Carlo and the DBFE method. The variance predicted by the DBFE method near the mean shock location is higher than the variance predicted by the Monte Carlo method for both velocity and pressure. As already observed for the Burgers equation, the higher variance in the DBFE method results from the oscillations near discontinuities associated with the Gibbs phenomenon. Effect of the Gibbs phenomenon is mitigated using the proposed postprocessing method, as can be observed from Fig. 9 through the close match in the variance between the Monte Carlo and the DBFE method with postprocessing.

Conclusion

This paper has 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 has been investigated for solution of a one-dimensional stochastic Burgers equation and stochastic quasi-one-dimensional flow through a convergent-divergent nozzle. The numerical results presented in this paper have demonstrated the ability of the proposed postprocessing method to mitigate the effect of the Gibbs phenomenon as discontinuities develop in the solution. Note that the proposed method does not require a priori 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.

Acknowledgment

This work was supported by AFOSR Grant (FA9550-12-1-0313).

2

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.

References

References
1.
Oreskes
,
N.
,
Shrader-Frechett
,
K.
, and
Belitz
,
K.
,
1994
, “
Verification, Validation and Confirmation of Numerical Models in Earth Sciences
,”
Science
,
263
(
5147
), pp.
641
647
.
2.
Mehta
,
U.
,
1991
, “
Some Aspects of Uncertainty in Computational Fluid Dynamics Results
,”
ASME J. Fluids Eng.
,
113
(
4
), pp.
538
543
.
3.
Mehta
,
U.
,
1996
, “
Guide to Credible Computer Simulations of Fluid Flows
,”
J. Propul. Power
,
12
(
5
), pp.
940
948
.
4.
Oberkampf
,
W.
,
DeLand
,
S.
,
Rutherford
,
B.
,
Diegert
,
K.
, and
Alvin
,
K.
,
2002
, “
Error and Uncertainty in Modeling and Simulation
,”
Reliab. Eng. Syst. Saf.
,
75
(
3
), pp.
335
357
.
5.
Walters
,
R.
, and
Huyse
,
L.
,
2002
, “
Uncertainty Analysis for Fluid Mechanics With Applications
,” National Aeronautics and Space Administration, Hanover, MD,
Report No. NASA/CR-2002-211449
.
6.
Rubinstein
,
R.
,
1981
,
Simulation and the Monte Carlo Method
,
Wiley
,
New York
.
7.
Chopin
,
N.
,
2009
, “
Central Limit Theorem for Sequential Monte Carlo Methods and Its Application to Bayesian Inference
,”
Ann. Stat.
,
32
(
6
), pp.
2385
2411
.
8.
Adler
,
R.
, and
Taylor
,
J.
,
2007
,
Random Fields and Geometry
,
Springer
,
New York
.
9.
Venturi
,
D.
,
2011
, “
A Fully Symmetric Nonlinear Biorthogonal Decomposition Theory for Random Fields
,”
Phys. D
,
240
(4–5), pp.
415
425
.
10.
Ghanem
,
R.
, and
Spanos
,
P.
,
2003
,
Stochastic Finite Elements: A Spectral Approach
,
Dover Publications
,
Mineola, NY
.
11.
Xiu
,
D.
, and
Karniadakis
,
G.
,
2003
, “
Modeling Uncertainty in Flow Simulations Via Generalized Polynomial Chaos
,”
J. Comput. Phys.
,
187
(
1
), pp.
137
167
.
12.
Xiu
,
D.
, and
Karniadakis
,
E.
,
2002
, “
The Weiner–Askey Polynomial Chaos for Stochastic Differential Equations
,”
SIAM J. Sci. Comput.
,
24
(
2
), pp.
619
644
.
13.
Koekoek
,
R.
, and
Swarttourw
,
R.
,
1998
, “
The Askey Scheme of Hypergeometric Orthogonal Polynomials and Its q-Analogue
,” Delft University of Technology, Delft, Netherlands,
Technical Report No. 98-17
.
14.
Mathelin
,
L.
,
Hussaini
,
M.
, and
Zang
,
T.
,
2005
, “
Stochastic Approaches to Uncertainty Quantification in CFD Simulations
,”
Numer. Algorithms
,
38
(1), pp.
209
236
.
15.
Najm
,
H.
,
2009
, “
Uncertainty Quantification and Polynomial Chaos Techniques in Computational Fluid Dynamics
,”
Annu. Rev. Fluid Mech.
,
41
(
1
), pp.
35
52
.
16.
LeVeque
,
R.
,
1992
,
Numerical Methods for Conservation Laws
,
Birkhause Verlag
,
Basel, Switzerland
.
17.
Sargsyan
,
K.
,
Safta
,
C.
,
Debusschere
,
B.
, and
Najm
,
H.
,
2012
, “
Uncertainty Quantification Given Discontinuous Model Response and a Limited Number of Model Runs
,”
SIAM J. Sci. Comput.
,
34
(
1
), pp.
44
64
.
18.
Najm
,
H.
,
2011
, “
Uncertainty Quantification in Fluid Flow
,”
Turbulent Combustion Modeling
(Fluid Mechanics and its Applications), Vol.
95
,
Springer, Dordrecht
,
The Netherlands
, pp.
381
407
.
19.
Hewitt
,
E.
, and
Hewitt
,
R.
,
1979
, “
The Gibbs–Wilbraham Phenomenon: An Episode in Fourier Analysis
,”
Arch. Hist. Exact Sci.
,
21
(
2
), pp.
129
160
.
20.
Gottlieb
,
D.
, and
Shu
,
C.
,
1997
, “
On the Gibbs Phenomenon and Its Resolution
,”
SIAM Rev.
,
39
(
4
), pp.
644
668
.
21.
Gottlieb
,
S.
,
Jung
,
J.
, and
Kim
,
S.
,
2011
, “
A Review of David Gottliebs Work on the Resolution of the Gibbs Phenomenon
,”
Commun. Comput. Phys.
,
9
(
3
), pp.
497
519
.
22.
Gottlieb
,
D.
,
Shu
,
C.
,
Solomonoff
,
A.
, and
Vandeven
,
H.
,
1992
, “
On the Gibbs Phenomenon—I: Recovering Exponential Accuracy From the Fourier Partial Sum of a Nonperiodic Analytic Function
,”
J. Comput. Appl. Math.
,
43
(1–2), pp.
81
98
.
23.
Gottlieb
,
D.
, and
Shu
,
C.
,
1995
, “
On the Gibbs Phenomenon—IV: Recovering Exponential Accuracy in a Subinterval From a Gegenbauer Partial Sum of a Piecewiae Analytic Function
,”
Math. Comput.
,
64
, pp.
1081
1095
.
24.
Sapsis
,
T.
, and
Lermusiaux
,
P.
,
2009
, “
Dynamically Orthogonal Field Equations for Continuous Stochastic Dynamical Systems
,”
Phys. D
,
238
(23–24), pp.
2347
2360
.
25.
Tagade
,
P.
, and
Choi
,
H.-L.
,
2012
, “
An Efficient Bayesian Calibration Approach Using Dynamically Biorthogonal Field Equations
,”
ASME
Paper No. DETC2012-70584.
26.
Cheng
,
M.
,
Hou
,
T.
, and
Zhang
,
Z.
,
2013
, “
A Dynamically Bi-Orthogonal Method for Time-Dependent Stochastic Partial Differential Equations—I: Derivation and Algorithms
,”
J. Comput. Phys.
,
242
, pp.
843
868
.
27.
Choi
,
M.
,
Sapsis
,
T.
, and
Karniadakis
,
G.
,
2013
, “
A Convergence Study for SPDEs Using Combined Polynomial Chaos and Dynamically-Orthogonal Schemes
,”
J. Comput. Phys.
,
245
, pp.
281
301
.
28.
Wan
,
X.
, and
Karniadakis
,
G.
,
2005
, “
An Adaptive Multi-Element Generalized Polynomial Chaos Method for Stochastic Differential Equations
,”
J. Comput. Phys.
,
209
(
2
), pp.
617
642
.
29.
Lin
,
G.
,
Su
,
C.
, and
Karniadakis
,
G.
,
2006
, “
Predicting Shock Dynamics in the Presence of Uncertainties
,”
J. Comput. Phys.
,
217
(
1
), pp.
260
276
.
30.
Poette
,
G.
,
Despres
,
B.
, and
Lucor
,
D.
,
2009
, “
Uncertainty Quantification for Systems of Conservation Laws
,”
J. Comput. Phys.
,
228
(
7
), pp.
2443
2467
.
31.
Tryoen
,
J.
,
Le Matre
,
O.
,
Ndjinga
,
M.
, and
Ern
,
A.
,
2010
, “
Roe Solver With Entropy Corrector for Uncertain Hyperbolic Systems
,”
J. Comput. Phys.
,
235
(
2
), pp.
491
506
.
32.
Chantrasmi
,
T.
,
Doostan
,
A.
, and
Iaccarino
,
G.
,
2009
, “
Pade–Legendre Approximants for Uncertainty Analysis With Discontinuous Response Surfaces
,”
J. Comput. Phys.
,
228
(
19
), pp.
7159
7180
.
33.
Cameron
,
R.
, and
Martin
,
W.
,
1947
, “
The Orthogonal Development of Non-Linear Functionals in Series of Fourier–Hermite Functionals
,”
Ann. Math.
,
48
(
2
), pp.
385
392
.
34.
Chakrabarti
,
A.
, and
Martha
,
S.
,
2009
, “
Approximate Solutions of Fredholm Integral Equations of the Second Kind
,”
Appl. Math. Comput.
,
211
(
2
), pp.
459
466
.
35.
Huang
,
S.
,
Quek
,
S.
, and
Phoon
,
K.
,
2001
, “
Convergence Study of the Truncated Karhunen–Loeve Expansion for Simulation of Stochastic Processes
,”
Int. J. Numer. Methods Eng.
,
52
(
9
), pp.
1029
1043
.
36.
Kurganov
,
A.
, and
Tadmor
,
E.
,
2000
, “
New High Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations
,”
J. Comput. Phys.
,
160
(
1
), pp.
241
282
.
37.
Gelb
,
A.
, and
Tadmor
,
E.
,
1999
, “
Detection of Edges in Spectral Data
,”
Appl. Comput. Harmonic Anal.
,
7
(
1
), pp.
101
135
.
38.
Gelb
,
A.
, and
Tadmor
,
E.
,
2001
, “
Detection of Edges in Spectral Data—II: Nonlinear Enhancement
,”
SIAM J. Numer. Anal.
,
38
(
4
), pp.
1389
1408
.
39.
Yang
,
H.
, and
Przekwas
,
A.
,
1992
, “
A Comparative Study of Advanced Shock-Capturing Schemes Applied to Burgers Equation
,”
J. Comput. Phys.
,
102
(
1
), pp.
139
159
.
40.
Pettersson
,
P.
,
Iaccarino
,
G.
, and
Nordstrom
,
J.
,
2009
, “
Numerical Analysis of the Burgers Equation in the Presence of Uncertainty
,”
J. Comput. Phys.
,
228
(
22
), pp.
8394
8412
.
41.
Mathelin
,
L.
,
Hussaini
,
M.
,
Zang
,
T.
, and
Bataille
,
F.
,
2004
, “
Uncertainty Propagation for a Turbulent, Compressible Nozzle Flow Using Stochastic Methods
,”
AIAA J.
,
42
(
8
), pp.
1669
1676
.