AN EFFICIENT COMPUTATIONAL METHOD FOR STATISTICAL MOMENTS OF BURGER’S EQUATION WITH RANDOM INITIAL CONDITIONS

The paper is concerned with efficient computation of numerical solutions to Burger’s equation with random initial conditions. When the Lax-Wendroff scheme (LW) is expanded using the Wiener chaos expansion (WCE), random and deterministic effects can be separated and we obtain a system of deterministic equations with respect to HermiteFourier coefficients. One important property of the system is that all the statistical moments of the solution to the Burger’s equation can be computed using the solution of the system only. Thus LW with WCE presents an alternative to computing moments by the Monte Carlo method (MC). It has been numerically demonstrated that LW with WCE approach is equally accurate but substantially faster than MC at least for certain classes of initial conditions.


Introduction
Uncertainty is observed in many and various phenomena in engineering, physics, biology, or finance.For example, small-scale effects in multiphase flow [13,14] or Navier-Stokes equations [18,26] may not be completely known, but subject to some random environmental effects.Or observational errors from inaccurate measurements may not be removed.Long waves of relatively shallow depth driven by white noise can be described approximately by Korteweg-de Vries equation with stochastic perturbations [8,30].Option pricing in financial markets may also introduce randomness [4,10].These uncertainties are represented mathematically by (functions of) random variables or stochastic processes.Then, the governing equation is given in the form of a stochastic partial differential equation (SPDE).Since the solution of this type of equation is random, one particular solution corresponding to a specific realization of a random variable is not of 2 Efficient computation of stochastic Burger's equation concern.Instead, it is important to know the statistical properties of the solution such as its mean, variance, or higher-order moments.
Wiener chaos expansion (WCE) is a Fourier expansion with respect to the randomness [3,15,19,23,24].In this paper, it is expanded with respect to the Gaussian random variable.The WCE of a solution to the Burger's equation, written by v(x,t,ξ), is given by where ξ n (x) is the nth-order Wick polynomial.As explained in Section 2, Wick polynomials are complete orthogonal basis in the Wiener space and the series in the right-hand side of (1.1) converges in L 2 (R,μ), where the Gaussian measure μ(dx) = (1/ √ 2π)exp{−x 2 / 2}dx.This polynomial chaos method has been used by many authors [5-7, 12, 17, 33].As every Fourier expansion will do, WCE separates variables.More specifically, it separates nonrandom variables (x,t) from the random one ξ.Then nonrandom Hermite-Fourier coefficients satisfy an infinite system of deterministic equations with nonlinearity similar to the one in Burger's equation.This system for the Hermite-Fourier coefficients is usually referred to as propagator, because it governs the propagation of randomness by the deterministic dynamics of the equation.One important property of the propagator is that all statistical moments of the solution to the Burger's equation can be computed by simple formulae that involve only the solution of the propagator.In addition, since the propagator is deterministic, it needs to be solved only once.Therefore, WCE presents an alternative to computing moments by the Monte Carlo (MC) method.There are two main issues addressed by the paper.First, a version of Lax-Wendroff scheme for the propagator of the Burger's equation has been developed.Second, the effectiveness of the WCE approach to computing statistical moments has been tested numerically and compared to the MC.A strong evidence has been demonstrated that suggested scheme is equally accurate but substantially faster than MC at least for certain initial conditions.
The Navier-Stokes (NS) equation is one of fundamental problems in fluid mechanics [20].When the turbulence, for example, is modeled using a random variable, the NS is expressed as a stochastic equation [15,18].Bensoussan and Temam initiated analytical study of an NS equation driven by a random force in [1].Then, [2,9,25,26] and many other researchers extended analytical studies of a stochastic Navier-Stokes (SNS) equation.But, there has not been many researches on numerical methods to solve an SNS.The MC is a classical way to solve an SPDE numerically [27,31].The MC may require a huge amount of computational time in order to obtain reliable estimation of statistical properties.Importance sampling or variance reduction method may be applied to reduce the workload of the MC [11,22].Selection of random number generators is another weakness of the MC.As a preliminary step to an SNS, an inviscid Burger's equation with a random initial condition is considered in the present study.As seen below, the method can be generalized to an SNS.See also [16,26].
The paper is organized as follows.Section 2 introduces Hermite polynomials and their properties.A version of Lax-Wendroff scheme for the propagator to the Burger's equation with random initial condition is derived in Section 3 using the WCE.Numerical results of the WCE approach to computing statistical moments are shown in Section 4.

Theoretical background
The inner product of two real-valued functions F(x) and G(x) is defined by F,G w ≡ ∞ −∞ F(x)G(x)w(x)dx, where the Gaussian probability density function, w(x) = (1/ √ 2π)e −x 2 /2 , is given as a weight function.Let G(x) w ≡ G,G w be the induced norm from the inner product.The Gram-Schmidt process with respect to the norm applied to the sequence of functions defined by 1,x,x 2 ,... (2.1) gives the orthonormal sequence {ξ n }, where (2.2) H n and ξ n are called the Hermite and Wick polynomials of order n, respectively.From the definition, the following relation can be easily derived: Theorem 2.1 (see [3,15]).{ξ n (x)} is a complete orthonormal set in L 2 (−∞,∞) with respect to the weight function w Reference [29] shows that a Fourier expansion in (1.1) converges uniformly.
Theorem 2.2 (see [29]).Let f (x) be defined for all finite values of x and integrable in any finite interval.Let F(x) = F(0) Since {H n } is a basis, the product H α (x)H β (x) can be represented with respect to {H n }, and the following lemma can be proved by induction.
where u α and v α are constants for all α.

Derivation of numerical algorithm
Let us consider the following conservation law for v(x,t,ξ) with a random initial condition: where f (v) = (1/2)v 2 and g(x,ξ) is a function of x and ξ. ξ is a Gaussian random variable N(0,1) with zero mean and unit variance.The WCE is a Fourier expansion with respect to the randomness.In this paper, it represents the solution v = v(x,t,ξ) of (3.1) as a series in ξ with respect to the Wick polynomials where v n is called the Hermite-Fourier coefficient of order n of v. H n (ξ) and ξ n (ξ) will be hereafter denoted by H n and ξ n , respectively, for notational simplicity.
When v(x,t,ξ) in (3.1) is expanded using a series (3.2), we obtain the following from Lemmas 2.3 and 2.4: Hongjoong Kim 5 Then, from the orthogonality in Theorem 2.1, for every n < ∞, the Hermite-Fourier coefficient v n is a solution of the following system: where t ∈ (0,T] and x ∈ R 1 .System (3.4) is referred to as propagator system because it governs the propagation of randomness.It separates nonrandom variables (x,t) from the random one ξ.Since this propagator system (3.4) is deterministic, the system needs to be solved only once.Once the propagator is solved, statistical moments can be computed as follows.First, note that E[ξ 0 ] = 1 and E[ξ n ] = 0 for n > 0, where E is the expectation symbol.Thus, the first moment is obtained by where the expectation and the limit can be interchanged because of the uniform convergence in Theorem 2.2.By Theorem 2.1 and Parseval's identity, the expression for the second moment can be derived as follows: In a similar way, higher moments can be obtained as functions of Hermite-Fourier coefficients.The Lax-Wendroff scheme (LW) for the deterministic conservation laws, where F j,i ≡ f (v j,i ) and A j,i ≡ (∂ f /∂v)(v j,i ) [28].Δ +x F j,i and Δ −x F j,i are defined by (F j+1,i − F j,i ) and (F j,i − F j−1,i ), respectively.A j±1/2,i denotes (1/2)(A j,i + A j±1,i ).When v j,i in (3.8) is expanded using a series in (3.2), the propagator system for each Hermite-Fourier coefficient v n j,i ≡ v n (x j ,t i ) can be solved by a following modified version of the Lax-Wendroff 6 Efficient computation of stochastic Burger's equation scheme (denoted by LW with WCE hereafter): Δt Δx where (3.10) , and 0 otherwise, and [x] is the smallest integer, which is not smaller than x.It is important to notice that (3.9) separates {z n,α j,i } and {w α j,i }, which removes unnecessary loops in computation.In addition, {w α j,i } get rid of repetitive calculations.Due to numerical limitations, algorithms with only a finite number of {v n } will be implemented.The number of {v n } used in the computation is called the length of the WCE.Equation (3.9) can be derived as follows: using (3.3), where Then, where where Hongjoong Kim 7 Then, (3.17) When (α,β,γ) are replaced by (p + q, p + n − q, p), (3.18) When (3.8) is expanded using these series, Δt Δx Δt Δx 8 Efficient computation of stochastic Burger's equation By comparing coefficients of ξ n , we obtain from orthogonality Δt Δx Δt Δx Since (3.20) contains nested loops, it is not practical to use this equation directly to update v n j,i .Instead, let α = p + n − q.Then, Note that the first term in the summation over α in the right-hand side of (3.21) is independent of m and k and that the second term is independent of q.This factorization removes unnecessary nested loops and saves computation.In addition, Notice also that since w α j,i is independent of n, it can be computed even prior to n-loop.Now by definition of z n,α j,i , completes the derivation of (3.9).

Hongjoong Kim 9
Initial and boundary conditions can be changed for v n as well.An initial condition v(x,0,ξ) = g(x,ξ) can be written as where g n = √ n!E[gξ n ] so that g n (x) becomes the initial condition for v n (x,t).Similarly, a boundary condition v(x 0 ,t,ξ) = h(t,ξ) can be written as where h n = √ n!E[hξ n ] and h n (t) is the boundary condition for v n (x,t) at x = x 0 .As [21] points out, for general systems of equations with arbitrary initial data, no numerical method has been proved to be stable or convergent in general, although convergence results have been obtained in some special cases.The stability of the LW with WCE in (3.9) is obtained by analyzing the stability of the scheme when the Jacobian of (3.1) is assumed to be a constant.Each Hermite-Fourier coefficient v n , then, satisfies deterministic Lax-Wendroff in (3.8) where a is the constant Jacobian so that the stability condition is where ν = (Δt/Δx)maxv is the Courant number.More analytical study on stability and error analysis will follow in the future research.As can be seen in the derivation, (3.9) can be generalized to the propagator of a viscous Burger's equation or even an SNS.

Numerical results
The stability condition (3.28) for the LW with WCE is numerically examined in Section 4 first.Then, accuracy of the scheme is investigated with respect to the length of the WCE and the Courant number.The accuracy is also compared to that from the MC.Finally, the analysis on efficiency of the scheme is shown based on the computation time.
where x ∈ [−1,1] and ξ ∼ N(0,1) is a Gaussian random variable with mean zero and variance one.For each realization of ξ, shocks will be developed.The initial condition can be obtained from (3.25): for x < 0, and for x > 0, Figure 4.1 shows the first and second moments of solutions for (4.1) computed by the LW with WCE when the Courant number is 0.6 or 1.0.Both solutions are stable, but the one obtained with 0.6 is much less accurate due to the known deficiency of the LW near discontinuities.If moments of v(x,t,ξ) are computed up to t ≤ T, as the Courant number is reduced, more computations will be required to reach the time T and the quality of the solution will be degraded, similarly to the deterministic LW [32].
It has been also found numerically that the stability condition becomes more conservative as the length of the WCE decreases, in the sense that the range of the Courant number for stability decreases as the length decreases.Table 4.1 shows maximum Courant numbers maintaining stability of the scheme for different lengths of the WCE.

Accuracy.
When (4.1) is solved, accuracy analysis is limited due to dispersions.Thus, let us consider the Burger's equation with the following random initial condition: where x ∈ [0,1] and ξ ∼ N(0,1) is Gaussian.From (3.25), Due to the nonlinearity, the statistical moments cannot be obtained analytically.When the MC is applied to the problem, statistical moments show fluctuations as the number of samples varies.The convergence has been obtained when the number of realizations exceeds 10 6 .Thus, the results from the MC with 10 6 samples will be considered as exact 12 Efficient computation of stochastic Burger's equation  values in the present study.For the error comparison, a relative error is used, which is defined by where u(x,t) is a statistical moment from the exact solution and u(x,t) is an approximation.Figure 4.2 shows relative errors of the MC in computing statistical moments of up to order 4 as the number of samples in the Monte Carlo simulation increases.Fluctuations of errors are observed in the figure even when the number of samples is large.fourth moments, the convergence is obtained with length 6, which is still small enough. Figure 4.8 shows relative errors from the LW with WCE for various Courant numbers.

Efficiency.
In order to evaluate the efficiency, let us compare computation times of these numerical schemes for the same relative error.For example, when the length of the WCE is 6 and the Courant number ν = 1.0, the LW with WCE results in 0.06% of relative error as seen in Figure 4.7.In order to obtain 0.06% of accuracy, Figure 4.2 shows that the MC requires more than 50000 realizations.the computational time, compared to the MC. Figure 4.10 also shows that as the Courant number decreases, Δt decreases and more computations are required.Another advantage of the LW with WCE over the MC is that once the propagator system is solved, since the output is a set of Hermite-Fourier coefficients, a statistical moment of any order can be estimated without any new computation.In case of the MC, however, orders of statistical moments of interest should be decided prior to the computation.For example, if first three moments were obtained during the MC computation, and if the fourth moment is needed after the computation is over, all the MC simulation should be performed again, because the fourth moment in this case cannot be generated from first three moments.

Conclusions
A version of LW has been derived using the WCE in order to compute efficiently statistical moments of solutions to Burger's equation with random initial conditions.The WCE transforms the LW into an infinite system of deterministic equations with respect to Hermite-Fourier coefficients.One of the important properties of the system is that all the statistical moments of the solution to the Burger's equation can be estimated by simple formulae that involve the solution of the system only.The stability, accuracy, and efficiency of the scheme have been numerically examined and compared to the MC.It has been demonstrated that LW with WCE is equally accurate but substantially faster than MC at least for certain classes of initial conditions.When the LW with WCE is applied to test problems, the computation time has been reduced up to 90%.In addition, the new approach does not require generation of random numbers and provides a good control over the computational errors.More study on stability and convergence will follow in the future research.

Figure 4 . 2 .
Figure 4.2.Errors of the MC in statistical moments as the number of simulations increases.

Figures 4 .
3 and 4.4 compare statistical moments of up to order 4 between the exact solution and the LW with WCE and show that very good accuracy can be obtained with sufficiently small lengths.Note also that dispersion or diffusion is not developed in Figures 4.3 and 4.4.Figures 4.5 and 4.6 show the differences between the statistical moments from the exact solution and those from the LW with WCE.Note that the scales in Figures 4.

5
and 4.6 are much finer than those in Figures 4.3

Figure 4 . 3 .
Figure 4.3.Comparison of the first (a) and second (b) moments of the exact solution with approximations from the LW with WCE for several lengths.

Figure 4 .Figure 4 . 4 .
Figure 4.7 shows relative errors (4.6) from the LW with WCE in statistical moments with respect to lengths of the WCE.As the length increases, more Hermite-Fourier coefficients are introduced and better accuracy is obtained.The first moment depends on v 0 only and good accuracy is obtained even when the length is 3.For the second, third, and

Figure 4 . 5 .
Figure 4.5.Differences in the first (a) and second (b) moments between the exact solution and approximations from the LW with WCE for several lengths.

Figure 4 .Figure 4 . 6 .Figure 4 . 7 .
Figure 4.6.Differences in the third (a) and fourth (b) moments between the exact solution and approximations from the LW with WCE for several lengths.

Figure 4 . 8 .
Figure 4.8.Relative errors of the LW with WCE in statistical moments for different Courant numbers.

Figure 4 .Figure 4 . 9 .Figure 4 . 10 .
Figure 4.10 shows the CPU time of the LW with WCE with respect to the lengths of the WCE when the Courant number is 0.6, 0.8, and 1.It takes 0.67 second with the LW with WCE when the length is 6 and the Courant number is 1.0, which is only 7.49% of

Table 4 .
1. Maximum Courant numbers maintaining stability of the LW with WCE.

Table 4 .
2 summarizes the efficiency of the LW with WCE over the MC.Relative error in the table is the error from the LW with WCE when the WCE of a given length is applied.

Table 4 .
2. Ratio of CPU times between the MC and the LW with WCE.MC and CPU LW in Table 4.2 represent the CPU times taken by the MC and the LW with WCE, respectively, to guarantee these relative errors.The table shows that the LW with WCE saves up to 90% of computational time in the given test problem. CPU