Stochastic Galerkin Projection and
Nu-merical Integration for Stochastic
Inves-tigations of the Viscous Burgers’
Equa-tion
Markus Wahlsten & Jan Nordstr¨om
Department of Mathematics
Link¨oping University
(will be inserted by the editor)
Stochastic Galerkin Projection and Numerical
Integration for Stochastic Investigations of the
Viscous Burgers’ Equation
Markus Wahlsten · Jan Nordstr¨om
Received: date / Accepted: date
Abstract We consider a stochastic analysis of the non-linear viscous Burg-ers’ equation and focus on the comparison between intrusive and non-intrusive uncer- tainty quantification methods. The intrusive approach uses a combina-tion of polynomial chaos and stochastic Galerkin projeccombina-tion. The non-intrusive method uses numerical integration by combining quadrature rules and the probability density functions of the prescribed uncertainties. The two meth-ods are applied to a provably stable formulation of the viscous Burgers’ equa-tion, and compared. As measures of comparison: variance size, computational efficiency and accuracy are used.
Keywords Uncertainty quantification· stochastic data · polynomial chaos · stochastic Galerkin· intrusive methods · non-intrusive methods · Burgers’ equation.
Mathematics Subject Classification (2000) 65D30· 65M06 · 35R60 · 35Q53
1 Introduction
The two main approaches for solving partial differential equations with random inputs are categorized in intrusive and non-intrusive methods. Semi-intrusive methods do exist, combining the intrusive and non-intrusive methods [1], [2], but are rare. Non-intrusive methods solves the original problem multiple times
Markus Wahlsten
Department of Mathematics, Computational Mathematics, Link¨oping University, SE-581 83 Link¨oping, Sweden (markus.wahlsten@liu.se).
Jan Nordstr¨om
Department of Mathematics, Computational Mathematics, Link¨oping University, SE-581 83 Link¨oping, Sweden (jan.nordstrom@liu.se).
using fix stochastic inputs [7], [5]. Numerical integration (NI) and interpola-tion techniques are then used to compute statistics of the soluinterpola-tion. Intrusive methods based on polynomial chaos (PC) expansions results in a system of equations for the expansion coefficients [14], [16]. The non-intrusive method unlike the intrusive one, typically relies on an already existing deterministic solver.
The aim of this article is to compare the performance of numerical NI and PC with the stochastic Galerkin approach. In general, a comparison between these two techniques are of course problem-dependent. However, focusing on problems related to fluid dynamics, we choose the viscous Burgers’ equation as a relevant model problem for the Navier–Stokes equations and compare the results in terms of computational time, variance size and accuracy.
The rest of the paper proceeds as follows. Section 2 introduces the non-linear viscous Burgers’ equation. Next, the fundamentals of polynomial chaos expansion is presented and the system of expansion coefficients is constructed and applied to the continuous problem, in Section 3. The foundation of (NI) is introduced in Section 4. In Section 5, we derive boundary conditions and an energy estimate for the continuous problem. Section 6 presents a stable and accurate semi-discrete finite difference formulation based on summation-by-parts (SBP) with simultaneous approximation terms (SAT). In Section 7, numerical results are presented comparing NI and PC. Finally, conclusions are drawn in Section 8.
2 The continuous problem
We consider the viscous Burgers’ equation in one space dimension ut+ uux = uxx, x∈ Ω, t > 0,
Lu = g(x, t, ξ), x∈ ∂Ω, t > 0, u = f (x, ξ), x∈ Ω, t = 0.
(1)
The solution is denoted u = u(x, t, ξ), where, ξ = (ξ1, ξ2, . . . , ξL) is the vector
of variables representing the uncertainty in the solution. The viscosity is a positive constant. The boundary operator defined on the boundary ∂Ω is de-noted by L, while f (x, ξ) and g(x, t, ξ) are the stochastic initial and boundary data.
3 Polynomial chaos expansion with stochastic Galerkin projection The foundation of polynomial chaos was first introduced in [6] and later gen-eralized in [15].
3.1 Polynomial chaos
Consider a general orthogonal stochastic basis{ψi(ξ)}∞i=0, satisfying
hψi, ψji =
Z
Ωξ
ψi(ξ)ψj(ξ) dP(ξ) = δij. (2)
whereP(ξ) is the probability density function of ξ. A second order random field u(x, t, ξ) satisfyingRΩξu(x, t, ξ)2dP(ξ) < ∞, can be expressed as
u(x, t, ξ) =
∞
X
k=0
uk(x, t)ψk(ξ). (3)
The coefficients{uk(x, t)}∞k=0 in (3) are given by the projections
uk(x, t) =hu(x, t, ξ), ψk(ξ)i, k = 0, 1, . . . (4)
By using (2)-(4), the expected value and variance can be expressed in terms of the polynomial chaos coefficients as
E[u] = Z
Ωξ
u dP(ξ) = u0, (5)
and
V ar[u] =E[u2]− E[u]2=
Z Ωξ ∞ X k=0 u2 kψk2dP(ξ) − u20= ∞ X k=1 u2 k. (6)
For more details on the polynomial chaos technique, see [10].
3.2 The stochastic Galerkin projection
Next, to perform a stochastic Galerkin projection of (1) we insert the truncated expansion u(x, t, ξ) = M X i=0 ui(x, t)ψi(ξ), (7) into (1), to get M X i=0 (ui)tψi+ M X i=0 M X j=0 ui(uj)xψiψj = M X i=0 (ui)xxψi, M X i=0 Luiψi= g(x, t, ξ), M X i=0 uiψi= f (x, ξ). (8)
By multiplying (8) with ψkand integrating over the stochastic domain Ωξ, we obtain (uk)t+ M X i=0 M X j=0 ui(uj)xhψiψj, ψki = (uk)xx, Luk=hg(x, t, ξ), ψki, uk=hf(x, ξ), ψki (9)
for k = 0, 1, . . . , M . Hence, a deterministic system of dimension M + 1 times the original system is obtained. From (9), the deterministic coefficients u0(x, t),
u1(x, t), . . . , uM(x, t) are computed.
4 Numerical integration
Numerical integration and quadrature relations in one dimension can be writ-ten Z b a f (q)ρ(q) dq≈ M X m=1 f (qm)wm, (10)
where f is the function to be integrated. The function ρ is the corresponding density function, the number of grid points is denoted M , and the quadrature points and weights are denoted qmand wmrespectively. For simplicity, we will
use the 4th-order accurate Simpson’s rule [4] as the integration technique in the rest of the paper.
A straightforward extension to higher dimensions is Z Γ f (˜q)ρ(˜q) d˜q≈ M1 X m1=1 · · · Mp X mp=1 f (qm1 1 , . . . , qpmp)wm1· · · wmp, (11)
where Γ is the p-dimensional domain and ˜q = (q1, . . . , qp). The quadrature
points and weights for dimension i are denoted qmi and wmi respectively. It is not necessary to use the same quadrature rule in every dimension. For high-dimensional problems adaptive sparse grid techniques [3] can be used to improve efficiency, however they are not considered in this paper.
5 The energy estimate
We now consider the system (9) written on the general form ut+ Aux= uxx, x∈ Ω, t > 0,
Lu = g, x∈ ∂Ω, t > 0, u = f, x∈ Ω, t = 0,
where we have chosen basis functions such that that A is symmetric and u = [u0, . . . , uM]T, Ajk= M X i=0 uihψiψj, ψki, g = [g0, . . . , gN], f = [f0, . . . , fN], gk=hg(x, ξ), ψki, fk=hf(x, ξ), ψki. (13)
The continuous problem (12) can be written in the split form [8] ut+
1
3(Au)x+
1
3Aux= uxx. (14)
By following the path in [11], we apply the energy method to obtain, kuk2t + 2kuxk2=− 2 3uTAu− 2uTux 1 0 =− u ux T2 3A−I −I 0 u ux 1 0, (15)
wherekuk2 =R01u2dx. Further, since the matrix A is symmetric, it can be
diagonalized in the following way
A(u) = XΛXT, X = [ ¯X, X 0], Λ = ¯Λ 0 0 Λ0 , ¯ X = [X+, X−], ¯Λ = Λ+ 0 0 Λ− . (16)
The matrices X+, X− and X0contains the positive, negative and zero
eigen-vectors of A respectively. Moreover, Λ+, Λ− and Λ0 are diagonal matrices
comprising the positive, negative and zero eigenvalues of A respectively. Remark 1 Due to non-linearity, the diagonalization in (16) depends on the solution, i.e. changes in time and space.
By diagonalizing (15), using (16) we get kuk2t+ 2kuxk2=−(W )TΛMW 1 0= BT, (17) where W = V u and V = ¯ XT− (2 3Λ)¯−1X¯T ∂∂x ¯XT ∂ ∂x I0− I0∂x∂ I0+ I0∂x∂ , ΛM= 2 3Λ¯ 0 0 0 0 −(2 3Λ)¯−1 0 0 0 0 I0/2 0 0 0 0 −I0/2 . (18) Note that in (18), potential zero eigenvalues are excluded since they do not contribute. The matrix I0is an identity matrix with the same size as Λ0. We
can further rewrite the right-hand side of (17) as BT =− W+ W− T Λ+M 0 0 Λ−M W+ W− x=1 + W+ W− T Λ+ M 0 0 Λ−M W+ W− x=0 (19)
where W+= V+u, W−= V−u, V+= X T +− (23Λ+)−1X+T∂x∂ XT −∂x∂ I0− I0∂x∂ , V−= X T −− (23Λ−)−1X−T∂x∂ XT +∂x∂ I0+ I0∂x∂ , Λ+M= 2 3Λ+ 0 0 0 −(2 3Λ−)−1 0 0 0 I0/2 , Λ− M= 2 3Λ− 0 0 0 −(2 3Λ+)−1 0 0 0 −I0/2 . (20)
The Cholesky decompositions of the positive definite matrices|Λ−M| and Λ+M are
Λ+M= R+RT+, |Λ−M| = R−RT−. (21)
The matrices R+,−in (21) are diagonal matrices containing the square root of
the eigenvalues of Λ+M and|Λ−M| respectively. By using (21) in (17) and (19) we get
kuk2t+ 2kuxk2=− (W+)T1(Λ+M)1(W+)1− (W−)T0|Λ−M|0(W−)0
+ (RT
−W−)T1(RT−W−)1+ (RT+W+)T0(RT+W+)0.
(22) To bound the left-hand side (LHS) of (22), we impose
(RT+W+)0= (RT+V+u)0= g0, (RT−W−)1= (R−TV−u)1= g1, (23)
where the subscript 0 and 1 denotes the left and right boundary, respectively. From (23) we can extract the boundary operators L0,1as
L0= (RT+V+)0, L1= (RT−V−)1. (24)
An example of this form of boundary conditions is described in A. By imposing the boundary conditions (23), (22) reduces to
kuk2t+ 2kuxk2=− (W+)T1(Λ+M)1(W+)1− (W−)T0|Λ−M|0(W−)0
+ gT
0g0+ g1Tg1. (25)
The result is summarized in
Proposition 1 The problem (12) augmented with the boundary conditions (23) has a bounded solution.
Proof By integrating (25) in time we obtain ku(T )k2+ 2 Z T 0 ku xk2dt≤ Z T 0 gT0g0+ gT1g1dt +kfk2. (26)
Hence, an energy estimate is obtained.
Remark 2 In NI, the above analysis is analogous, the only differences is that u = u(x, t, ξ), g0,1= g0,1(t, ξ), f = f (x, ξ) and A = u.
6 The semi-discrete formulation
The problem (9) or equivalently (12) is solved using a finite difference formu-lation based on the SBP–SAT technique [12], [13], [9]. A stable and accurate semi-discrete formulation of (12) on SBP–SAT form using the split form (14) is vt+13(D⊗ IM)Av + 31A(D⊗ IM)v− (D2⊗ IM)v + (P−1E 0⊗ IM)Σ0(L0v− e0⊗ ˜g0) + (P−1EN⊗ IM)Σ1(L1v− eN x⊗ ˜g1) v(0) = f , (27)
where⊗ denotes the Kronecker product, ˜g0= [g0T, 0]T and ˜g1= [0, gT1]T. In
(27), v represents u numerically, g0, g1and f are the numerical approximations
ofhg0(t, ξ), ψli, hg1(t, ξ), ψli, and hf(x, ξ), ψli for l = 0, 1, . . . , M, respectively.
The matrices Σ0and Σ1are to be determined to ensure stability.
The numerical solution is arranged in the following way
v = v0 v1 .. . vi .. . vNx , vi= vi0 vi1 .. . vim .. . viM ,
where vimapproximates the polynomial chaos coefficient um(xi, t).
Remark 3 The vector vim approximates u(xi, t, ξm) when computing the
so-lution using NI.
The approximate derivative in the x−direction is approximated by the SBP operator D = P−1Q. The matrix P is a positive definite diagonal
ma-trix and Q is almost skew-symmetric satisfying Q + QT = E
N− E0 =B =
diag[−1, 0, . . . , 0, 1]. We denote the identity matrix of dimension M + 1 by IM. The matrices E0 and EN are zero except for the element (1, 1) and
(Nx+ 1, Nx+ 1) respectively, which is 1. Finally, eNx denotes a zero vector with the exception of the last element which is 1.
The boundary and initial data g0, g1and f , consist of the projections
g0,1= h¯g0,1(t, ξ), ψ0(ξ)i .. . h¯g0,1(t, ξ), ψM(ξ)i , f = h ¯f (x0, ξ), ψ0(ξ)i .. . h ¯f (x0, ξ), ψM(ξ)i .. . h ¯f (xNx, ξ), ψ0(ξ)i .. . h ¯f (xNx, ξ), ψM(ξ)i . (28)
The original boundary and initial data vectors as a function of ξ are ¯gi(t, ξ)
and ¯f (x, ξ) respectively, injected on the boundaries x = i and t = 0. The inner productsh¯g0(t, ξ), ψm(ξ)i, h¯g1(t, ξ), ψm(ξ)i and h ¯f (x, ξ), ψm(ξ)i are computed
numerically using NI. The matrix A in (27) is given by A = diag( ¯A0, . . . , ¯ANx), ( ¯Ai)jk=
M
X
l=0
vilhψl, ψj, ψki. (29)
Remark 4 When using NI, the vectors g0, g1 and f instead denote g0,1 =
[¯g0,1(t, ξ0), . . . , ¯g0,1(t, ξM)]Tand f = [ ¯f (x0, ξ0), . . . , ¯f (x0, ξM), . . . , ¯f (xNx, ξ0), . . . , ¯
f (xNx, ξM)]T. Hence, the boundary and initial data ¯g0(ξ), ¯g1(ξ) and ¯f (ξ) are grid functions in ξ. Moreover, the matrix blocks ¯Ai in the NI framework
cor-respond to
( ¯Ai) = diag(vi0, . . . , viM). (30)
The discrete boundary operators L0and L1in (24) are decomposed as
L0= [RT+V+], L1= [RT−V−]. (31) where V−,+0 = X T −,+ 0 I0 , V−,+D = − 3 2(Λ−,+)−1XT−,+ XT −,+ I0 , V−,+= V−,+0 + V−,+D D, Λ˜ −,+M = 2 3Λ−,+ 0 0 0 −3 2(Λ+,−)−1 0 0 0 ∓I0/2 , (32)
and ˜D = (D⊗ IM). Further, Λ+M and|Λ−M| are positive definite and can be
split as Λ+M= R+RT+= q Λ+M q Λ+M, |Λ− M| = R−RT−= q |Λ− M| q |Λ− M| (33)
In (32), the matrices X+and X−contain the eigenvectors of their
correspond-ing positive and negative eigenvalues of A respectively.
6.1 Stability
The discrete energy method (multiplying (27) by vT(P⊗ I
M) and adding the
transpose), yields d dtkvk 2 P + 1 3v T((Q ⊗ IM)A + A(QT⊗ IM))v − 13vT(A(Q⊗ I M) + (QT⊗ IM)A)v = vT(Q⊗ I M) ˜Dv + ( ˜Dv)T(QT⊗ IM)v + 2vT(E 0⊗ IM)Σ0(L0v− e0⊗ ˜g0) + 2vT(E N⊗ IM)Σ1(L1v− eN⊗ ˜g1). (34)
The SBP-property applied to (34) together with the splits Σ0= diag(Σ+0, 0)
and Σ1= diag(0, Σ−1) results in
d dtkvk 2 P + 2 ˜Dv 2 P = v ˜Dv T 0 2 3A−I −I 0 0 v ˜Dv 0 − v ˜Dv T 1 2 3A−I −I 0 1 v ˜Dv 1 + 2((Σ+0)Tv)T 0((RT+V+v)0− g0) + 2((Σ−1)Tv)T 1((RT−V−v)1− g1). (35) where we denote kvk2P = vT(P ⊗ I M)v and ˜Dv 2 P = ( ˜Dv) T(P ⊗ I M) ˜Dv.
Further, we make the choices Σ+0 = (RT
+V+)TΣ0, Σ−1 = (RT−V−)TΣ1, (36)
and diagonalize A = XΛXT as in the continuous setting (16) to obtain
d dtkvk 2 P+ 2 ˜Dv 2 P =− V+v V−v T 1 Λ+M 0 0 Λ−M 1 V+v V−v 1 + V+v V−v T 0 Λ+M 0 0 Λ−M 0 V+v V−v 0 + RT +V+v g T 0 2Σ0 −Σ0 −Σ0 0 0 RT +V+v g 0 + RT −V−v g T 1 2Σ1 −Σ1 −Σ1 0 1 RT −V−v g 1 . (37)
We rewrite (37) and obtain d dtkvk 2 P+ 2 ˜Dv 2 P =− (W −)T 0(|Λ−M|)0(W−)0− (W+)T1(Λ+M)1(W+)1 + RT +W+ g T 0 I + 2Σ0−Σ0 −Σ0 0 0 RT +W+ g 0 + RT −W− g T 1 I + 2Σ1−Σ1 −Σ1 0 1 RT −W− g 1 , (38) where we have used the split Λ+M = R+RT+,|Λ−M| = R−RT−and introduced
the notations W−,+= V−,+v.
The choice
and adding and subtracting gT 0g0and g1Tg1we get d dtkvk 2 P+ 2 ˜Dv 2 P =− (W −)T 0(|Λ−M|)0(W−)0 − (W+)T 1(Λ + M)1(W+)1 − RT +W+ g T 0 (C⊗ I) RT +W+ g 0 − RT −W− g T 1 (C⊗ I) RT −W− g 1 + gT 0g0+ gT1g1 (40) where C = 1−1 −1 1 . (41)
The matrix C has the eigenvalues 0, 2, and is hence positive semi-definite, which in turn means that the LHS of (40) is bounded by data.
We summarize the results in
Proposition 2 The numerical approximation (27) with the penalty matrices Σ0= diag(−(RT+V+)T, 0), Σ1= diag(0,−(RT−V−)T) (42)
is strongly stable.
Proof By integrating (40) in time from 0 to T we obtain kv(T )k2P + 2 Z T 0 ˜Dv 2 P dt≤ Z T 0 g0Tg0+ gT1g1dt +kfk2P. (43)
The LHS of (43) is bounded by data, leading to strong stability. Note the re-semblance with the continuous energy estimate (26).
7 Numerical experiments
To exemplify the difference between NI and PC, we consider the split form (14) of (12) with boundary conditions of the form (24). The initial and boundary data are given by the manufactured solution
w(x, t, ξ) = 5 + e−ξ/2sin(µπe−ξ2/2x− t). (44)
In (44), an increased µ leads to an increased variation of the stochastic solution. The model problem considered for the numerical experiments is (1) aug-mented with the forcing function F = wt+ wwx− wxx, and = 0.01. As a
measure of comparison, the following quantity is used V arV =
Z T
0 kV ar[U] − V ar[V ]k2
In (45), U denotes a PC computation using 25 basis functions (on the same spatial grid) is used in order to reduce the deterministic errors. U is considered to be the ”exact” solution, while V denotes the computed numerical solution. The comparison is done using µ = 1 and µ = 5. A small µ is generates a slow varying solution. In the calculations below a 3rd-order SBP-operator with 40 grid points in space, and a 4th-order Runge–Kutta scheme in time is used.
Figure 1a and 1b illustrate the error of the variance as a function of num-ber of coefficients/evaluations (M ) for PC and NI using µ = 1 and µ = 5, respectively. Figure 1c and 1d show the CPU time as a function of M for PC and NI for the same two cases. Finally, Figure 1e and 1f depict the error of the variance as a function of CPU time for PC and NI using again µ = 1 and µ = 5.
As can be seen from Figure 1a-1f, PC performs better than NI for slow varying problems, while NI seems to be more efficient for fast varying problems.
8 Summary and conclusions
We have analyzed and compared the efficiency of PC and NI. The study has been carried out on the non-linear viscous Burgers’ equation, where non-linear boundary conditions were derived resulting in an energy bound and non-linear stability. The PC framework is employed to the continuous problem, and a sta-ble high-order finite difference formulation on SBP-SAT form was constructed. A similar numerical scheme was constructed for NI.
The difference in the variance was used as a measure of comparison. The numerical results suggest that the PC procedure outperforms NI for slow vary-ing problems, while NI seems to be more efficient for fast varyvary-ing problems. The difference in performance opens the door for possible gains using hybrid methods.
A A boundary operator example
To show a simple example of the form of the boundary operators, we let
A = u1u2 u2u1 (46)
which can be diagonalized as
A = XΛXT, X = √1 2 −1 1 1 1 , Λ = u1− u2 0 0 u1+ u2 . (47)
The matrices in (47) are decomposed into
Λ = Λ+ 0 0 Λ− , X = [X+, X−]. (48)
100 101 102 M 10-10 10-5 100 NI PC
(a) The error of the variance as a function of M using µ = 1. 100 101 102 M 10-8 10-6 10-4 10-2 100 NI PC
(b) The error of the variance as a function of M using µ = 5. 100 101 102 M 101 102 CPU time (s) NI PC
(c) The total CPU time as a function of M using µ = 1. 100 101 102 M 101 102 CPU time (s) NI PC
(d) The total CPU time as a function of M using µ = 5. 101 102 CPU time (s) 10-10 10-5 100 NI PC
(e) The error of the variance as a function of the CPU time using µ = 1.
101 102 CPU time (s) 10-8 10-6 10-4 10-2 100 NI PC
(f) The error of the variance as a function of the CPU time using µ = 5.
To avoid a complicated structure, we assume at the boundaries x = 0, 1, that u1< u2
and u1, u2> 0, which leads to
V+= √1 2 " 1− 3 2(u1+u2) ∂ ∂x 1− 3 2(u1+u2) ∂ ∂x −∂ ∂x ∂ ∂x # V−= √1 2 " −1 + 3 2(u1−u2) ∂ ∂x 1− 3 2(u1−u2) ∂ ∂x ∂x∂ ∂x∂ # (49) and R−= q 2 3|u1− u2| 0 0 q32u 1 1+u2 , R+= q 2 3(u1+ u2) 0 0 q3 2 1 |u1−u2| . (50) Finally, by combining (49) and (50) we obtain the boundary operators
L0= (RT+V+)0 = q u1+u2 3 − 2 q 3 u1+u2 ∂ ∂x q u1+u2 3 − 2 q 3 u1+u2 ∂ ∂x − 2 q 3 |u1−u2| ∂ ∂x 2 q 3 |u1−u2| ∂ ∂x 0 , L1= (RT−V−)1 = − q |u1−u2| 3 − 2 q 3 |u1−u2| ∂ ∂x q |u1−u2| 3 + 2 q 3 |u1−u2| ∂ ∂x 2 q 3 u1+u2 ∂ ∂x 2 q 3 u1+u2 ∂ ∂x 1 . (51) References
1. Abgrall, R., Congedo, P.M.: A semi-intrusive deterministic approach to uncertainty quantification in non-linear fluid flow problems. J COMPUT PHYS 235, 828–845 (2013)
2. Abgrall, R., Congedo, P.M., Corre, C., Galera, S.: A simple semi-intrusive method for uncertainty quantification of shocked flows, comparison with a non-intrusive polynomial chaos method. In: ECCOMAS CFD, pp. 1–11 (2010)
3. Bungartz, H.J., Dirnstorfer, S.: Multivariate quadrature on adaptive sparse grids. Com-puting 71(1), 89–114 (2003)
4. Davis, P.J., Rabinowitz, P.: Methods of numerical integration. Courier Corporation (2007)
5. Eldred, M., Burkardt, J.: Comparison of non-intrusive polynomial chaos and stochastic collocation methods for uncertainty quantification. AIAA paper 976(2009), 1–20 (2009) 6. Ghanem, R.G., Spanos, P.D.: Stochastic finite elements: a spectral approach. Courier
Corporation (2003)
7. Hosder, S., Walters, R.W., Perez, R.: A non-intrusive polynomial chaos method for uncertainty propagation in CFD simulations. AIAA paper 891, 2006 (2006)
8. Nordstr¨om, J.: Conservative finite difference formulations, variable coefficients, energy estimates and artificial dissipation. J SCI COMPUT 29(3), 375–404 (2006)
9. Nordstr¨om, J., Carpenter, M.H.: Boundary and interface conditions for high order finite difference methods applied to the Euler and Navier-Stokes equations. J COMPUT PHYS 148, 621–645 (1999)
10. Pettersson, M.P., Iaccarino, G., Nordstr¨om, J.: Polynomial chaos methods for hyperbolic partial differential equations: Numerical techniques for fluid dynamics problems in the presence of uncertainties. Springer (2015)
11. Pettersson, P., Iaccarino, G., Nordstr¨om, J.: Numerical analysis of the Burgers equation in the presence of uncertainty. Journal of Computational Physics 228(22), 8394–8412 (2009)
12. Strand, B.: Summation by parts for finite difference approximations for d/dx. J COM-PUT PHYS 110(1), 47–67 (1994)
13. Sv¨ard, M., Nordstr¨om, J.: On the order of accuracy for difference approximations of initial-boundary value problems. J COMPUT PHYS 218(1), 333–352 (2006) 14. Wan, X., Karniadakis, G.E.: An adaptive multi-element generalized polynomial chaos
method for stochastic differential equations. J COMPUT PHYS 209(2), 617–642 (2005) 15. Xiu, D., Karniadakis, G.E.: Modeling uncertainty in steady state diffusion problems via generalized polynomial chaos. COMPUT METHOD APPL M 191(43), 4927–4948 (2002)
16. Xiu, D., Karniadakis, G.E.: The Wiener–Askey polynomial chaos for stochastic differ-ential equations. SIAM J SCI COMPUT 24(2), 619–644 (2002)