A Study of The Stochastic Burgers’ Equation Using The Dynamical Orthogonal Method

Abstract

In the current work, the stochastic Burgers’ equation is studied using the Dynamically Orthogonal (DO) method. The DO presents a low-dimensional representation for the stochastic fields. Unlike many other methods, it has a time-dependent property on both the spatial basis and stochastic coefficients, with flexible representation especially in the strongly transient and nonstationary problems. We consider a computational study and application of the DO method and compare it with the Polynomial Chaos (PC) method. For comparison, both the stochastic viscous and inviscid Burgers’ equations are considered. A hybrid approach, combining the DO and PC is proposed in case of deterministic initial conditions to overcome the singularities that occur in the DO method. The results are verified with the stochastic collocation method. Overall, we observe that the DO method has a higher rate of convergence as the number of modes increases. The DO method is found to be more efficient than PC for the same level of accuracy, especially for the case of high-dimensional parametric spaces. The inviscid Burgers’ equation is analyzed to study the shock wave formation when using the DO after suitable handling of the convective term. The results show that the sinusoidal wave shape is distorted and sharpened as the time evolves till the shock wave occurs.

1. Introduction

Recently, there has been a compelling need to quantify uncertainties of parameters in physical and mathematical problems in a probabilistic framework. This includes cases that can be characterized by stochastic partial differential equations (SPDEs). Theses equations appear in many application fields, for example, solid mechanics, random vibrations, fluid dynamics [1,2,3], propagation of waves through random media [4,5,6], and finance [7]. Stochastic forcing, uncertainties in one or more physical parameters, and initial and/or boundary values, among others, are contributing to stochasticity in many disciplines.

Burgers’ equation is one of many important models that appear in many applications specially in fluid mechanics. When considering the stochastic effects added to Burgers’ equation, it will be extended to study some real-life applications such as the flow turbulence. There are various techniques that can be used to study and analyze such stochastic models. The Monte Carlo (MC) simulation [8,9,10] is one of the most practical methods, using a sequence of random numbers to handle these problems and get useful statistical properties. To attain a certain accuracy level, a sufficient number of samples are required. More studies have been applied to improve the MC efficiency such as sequential Monte Carlo methods [11,12] and the Error Subspace Statistical Estimation (ESSE) [13,14]. The MC-based techniques still have a relatively low rate of convergence and/or do not work efficiently for the general case [15].

Order-reduction techniques have been employed to simplify and analyze high-dimensional complicated systems relative to several scientific and engineering problems. Hence, they provide lower complexity for the original SPDE model. One of these approaches is the Proper Orthogonal Decomposition (POD) [16], which is based on the statistical technique of Karhunen-Loeve (KL) expansion. The main drawback is that the spatial basis of the POD method is selected a priori, which means that it may be unable to describe some problems, such as transient fluid flows, which are highly time-dependent and make the basis used irrelevant as time evolves.

Ghanem and Spanos [17] introduced an important and widely used approach known as the Polynomial Chaos (PC) expansion, which is based on Wiener’s theory on polynomial chaos [18]. The main drawback is that the PC suffers from problems relying on long time dependence, although it is well known that the PC method shows fast convergence rates for Gaussian processes. Other types of processes may have slower convergence rates in addition to solution deterioration with time evolution. The generalized PC expansion (gPC) introduced by Xiu and Karniadakis [19,20] uses different basis functions for different types of problems and is employed to increase the rate of convergence and to improve the efficiency of a wide range of nonlinear applications. The Probabilistic Collocation Method (PCM) [21] is an efficient version derived from PC that provides a smooth solution in parametric spaces and gives fast convergence rates as the order of the expansion increases [22,23].

The PC can also be extended to analyze the stochastic differential equations (SDEs) which are derived by or include white noise. The Wiener chaos expansion (WCE) and the Wiener-Hermite expansion (WHE) are common techniques used for analyzing SDEs [24]. The fractional-order version of the WHE was recently developed for the case of fractional Brownian motion [25].

In the current study, we shall focus on analyzing Burgers’ model with stochastic forcing. The stochastic Burgers’ model is common in studying flow turbulence and in general the propagation of randomness and nonlinear waves in the undispersive media [24]. The same type of nonlinearity of Burgers’ equation exists as in Navier-Stokes equations. Results obtained from solving the 1D stochastic Burgers’ equation are similar to experimental results of the 3D real turbulence problems.

Dynamical Orthogonal (DO) decomposition is introduced in [26,27] as one of the reduced-order techniques in which the solution is approximated in a generalized KL expansion, i.e.,

$$u\left(x,t;\omega \right)=\overline{u}\left(x,t\right)+{\displaystyle \sum }_{i=1}^N{u}_i\left(x,t\right)Y_i\left(t;\omega \right),$$

where x is an n-dimension space vector, t > 0 is the time, $\omega$ is the output of a random experiment, $N$ represents the number of approximating modes,$\overline{u}\left(x,t\right)$ is the mean field, and $u_i\left(x,t\right) ;i=0 to N,$ are orthonormal time-dependent fields in the spatial domain. The stochastic processes $Y_i\left(t;\omega \right)$ are usually considered to be with zero mean without loss of generality. Both spatial basis and the stochastic coefficients have a time dependency. As a result, the above KL representation is particularly flexible when it comes to represent nonstationary, highly transient responses. On the other hand, the same property adopts a redundant representation. To overcome this redundancy in both stochastic coefficients and spatial basis, a natural constraint is imposed, called the dynamical orthogonal condition [26]. By means of this condition, the DO components, i.e., $u_i\left(x,t\right),Y_i\left(t;\omega \right)$ and $\overline{u}\left(x,t\right)$, can be derived.

From a computational viewpoint, the uncertainty evolution using DO decomposition is accomplished by determining a set of $\left(N+1\right)$ deterministic PDEs that describe the mean field evolution $\overline{u}\left(x,t\right)$ and the basis $u_i\left(x,t\right)$ coupled with $N$ (ordinary) SDEs that describe the evolution of stochastic coefficients $Y_i\left(t;\omega \right)$. The DO equations reach a singularity at very low levels of uncertainty when the modes develop independently. Therefore, this may generate a numerical issue since the computations include ratios of small moment quantities, particularly in problems where deterministic initial conditions are included. The initialization of the random space is another issue that arises when starting with deterministic initial conditions.

To address these issues, a hybrid approach is introduced in [28] that avoids the singularity that occurs in the covariance matrix. This hybrid approach was formulated by integrating the PC with the DO method. Firstly, the PC is used in solving the SPDE up to a certain time until guaranteeing the development of stochasticity. Next, the DO approach is used to continue the solution with time evolution. As part of this strategy, the KL expansion is employed to give a collection of modes that are used for initializing the stochastic coefficients.

In the current work, we study the stochastic Burgers’ model using the DO method and compare it with the PC from different sides including error behavior, effect of switching time, and the number of modes and their effect on solution accuracy. Moreover, the error between DO components and the growth of eigenvalues are computed. The DO behavior for the case of inviscid Burgers’ equation and shock wave formation is also analyzed in the current work. Suitable known techniques, such as upwind techniques [29], are combined with the DO method to analyze the shock formation in the presence of random coefficients.

The paper is structured as follows: In Section 2, we provide a brief overview of the construction of the DO representation, the evolution equations, and their numerical solution. A hybrid approach combining the DO and PC is presented. In Section 3, the performance of the DO and PCM methods are compared by analyzing the stochastic viscous and inviscid Burgers’ equations. In Section 4, we study shock wave formation via inviscid Burgers’ equation. In Section 5, we present the conclusions.

2. Methodology

Consider the probability space ($\mathsf{\Omega } , \mathcal{F},P$) where the sample space $\mathsf{\Omega }$ contains all the elementary events $\omega . \mathcal{F}$ is the $\sigma$-algebra of the subsets of $\mathsf{\Omega }$, and $P$ is the probability measure. The random field will be defined for all measurable maps with the form $u\left(x,t,\omega \right)$. The mean value operator of $u$ is formulated as follows:

$$\overline{u}\left(x,t\right)=E\left[u\left(x,t;\omega \right)\right]={{\displaystyle \int }}_{\mathsf{\Omega }}u\left(x,t;\omega \right)dP\left(\omega \right).$$

(1)

The set of continuous random fields that are square integrable, i.e., ${{\displaystyle \int }}_{\mathrm{D}}E[u{\left(x,t;\omega \right)}^{T}u\left(x,t;\omega \right)] dx<\infty$, where the transpose of $u\left(x,t;\omega \right)$ is $u{\left(x,t;\omega \right)}^T$ for all $t ϵ T$; where D is the domain and T is the time interval. The covariance operator, between two random fields $u\left(x,t;\omega \right)$ and $q\left(x,t;\omega \right)$, has the bi-linear form:

$$C_{u, q}\left(x,y\right)=E\left[{\left(u\left(x,t;\omega \right)-\overline{u}\left(x,t\right)\right)}^T\left(q\left(y,\tau ;\omega \right)-\overline{q}\left(y,\tau \right)\right)\right],$$

(2)

where $x,y \mathrm{belong} \mathrm{to} \mathrm{the} \mathrm{spatial} \mathrm{domain} D, t,\tau \in T,$ and forms a Hilbert space H [4,30]. For the random $u\left(x,t;\omega \right)$ the projection operator $\prod$ to subspace spanned by the orthonormal basis $Q={\left\{s_j\left(x,t;\omega \right)\right\}}_{j=1}^m$ is defined as follows:

$${{\displaystyle \prod }}_Q\left[u\left(x,t;\omega \right)\right]={\displaystyle \sum }_{\mathrm{j}=1}^m{s}_j\left(.,t;\omega \right),u\left(.,t;\omega \right) s_j\left(x,t;\omega \right).$$

(3)

We define the autocovariance operator for $u\left(x,t;\omega \right)$ in the case $\tau =t$ as follows:

$$C_{u,u}\left(x,y\right)=E\left[{\left(u\left(x,t;\omega \right)-\overline{u}\left(x,t\right)\right)}^T\left(u\left(y,t;\omega \right)-\overline{u}\left(y,t\right)\right)\right],$$

(4)

where $x,y\in D, t\in T$. Then, we define the integral operator at a given time $t$ which depends on the autocovariance operator as follows:

$${\mathsf{{\rm I}}}_C \varphi ={{\displaystyle \int }}_D{C}_{u,u}\left(x,y\right) \varphi \left(x,t\right)dx\hspace{1em};\varphi \in L^2.$$

(5)

and it is a positive operator, compact and self-adjoint in the Hilbert space $H$ of a deterministic, square-integrable, and continuous field $L^2$ [16,31].

In the current work, we shall consider the following evolution SPDE:

$$\frac{\partial u\left(x,t;\omega \right)}{\partial t}=\mathcal{L}\left[u\left(x,t;\omega \right),\omega \right]\hspace{1em}; x\in D, t\in T, \omega \in \mathsf{\Omega },$$

(6)

where $\mathcal{L}$ is the differential operator, generally nonlinear. The system initial condition at time $t_0$ is given by the random field:

$$u\left(x,t_0; \omega \right)=u_0\left(x;\omega \right).$$

(7)

and the boundary conditions take the form:

$$B\left[u\left(\xi ,t;\omega \right)\right]=h\left(\xi ,t;\omega \right)\hspace{1em}; \xi \in \partial D,$$

(8)

where $B$ is a differential linear operator and $\partial D$ is the boundary of domain $D$ in $R^n$ and $n=1, 2$ or $3$.

2.1. DO Method

For every random field $u\left(x,t,\omega \right)\in L^2$, the KL expansion [29,30] at time $t$ can be written as:

$$u\left(x,t;\omega \right)=\overline{u}\left(x,t\right)+{\displaystyle \sum }_{i=1}^{\infty }{Y}_i\left(t;\omega \right) u_i\left(x,t\right),$$

(9)

where $u_i\left(x,t\right)$ are the eigenfunctions and $Y_i\left(t;\omega \right)$ are stochastic processes. We can obtain the eigenfunctions and eigenvalues of the covariance matrix $C$ of $u_i\left(x,t\right)$ by solving:

$${{\displaystyle \int }}_D{C}_{u, u}\left(x,y\right) u_i\left(x,t\right) dx={\lambda }_i u_i\left(y,t\right),$$

(10)

where $u_i$ are orthonormal eigenfunctions, i.e., $u_i,u_j={\delta }_{ij}$. Then, we can approximate the random field $u\left(x,t,\omega \right)$ by a finite N-modes series expansion as:

$$u\left(x,t,\omega \right)=\overline{u}\left(x,t\right)+{{\displaystyle \sum }}_{i=1}^N{Y}_i\left(t;\omega \right) u_i\left(x,t\right).$$

(11)

The stochastic subspace $V_s=span{\left\{u_i\left(x,t\right)\right\}}_{i=1}^N$ is defined as the space spanned by the eigenfunctions of the largest $N$ eigenvalues. It is worth noting that the spatial basis $u_i\left(x,t\right)$ and stochastic coefficients $Y_i\left(t;\omega \right)$ are recognized to have a time-dependent nature, and they change in relation to system dynamics.

The mean component $\overline{u}\left(x,t\right),$ the orthogonal basis ${\left\{u_i\left(x,t\right)\right\}}_{i=1}^N$, and the stochastic coefficients ${\left\{Y_i\left(t;\omega \right)\right\}}_{i=1}^N$ are all time-dependent and they are varyingly redundant. The main cause of redundancy arises from the uncertainty of evolution. In order to solve this redundancy, we can limit the basis ${\left\{u_i\left(x,t\right)\right\}}_{i=1}^N$ evolution to be only orthogonal to the space $V_s$ by imposing additional constraints. For example, we can apply the following orthogonality constraint:

$$\frac{d{V}_s}{dt} \perp V_s\iff \langle \frac{\partial u_j\left(x,t\right)}{\partial t}, u_i\left(x,t\right)\rangle =0\hspace{1em}\hspace{1em}; i,j=1,\dots ,N.$$

(12)

This constraint is known as the dynamical orthogonal condition [26,32]. This is also equivalent to preserving the basis ${\left\{u_i\left(x,t\right)\right\}}_{i=1}^N$ orthogonality for all times since:

$$\frac{\partial }{\partial t}\langle u_j\left(.,t\right),u_i\left(.,t\right)\rangle =\langle \frac{\partial u_j\left(x,t\right)}{\partial t}, u_i\left(x,t\right)+u_j\left(x,t\right),\frac{\partial u_i\left(x,t\right)}{\partial t}\rangle .$$

(13)

The DO expansion is used to obtain a set of independent equations characterizing all deterministic and stochastic quantities. Using DO, we will be able to reduce the SPDE to a set of $N+1$ deterministic PDEs for $\overline{u}\left(x,t\right)$ and ${\left\{u_i\left(x,t\right)\right\}}_{i=1}^N$ coupled with $N$ ordinary differential equations in the stochastic coefficients ${\left\{Y_i\left(x,t\right)\right\}}_{i=1}^N$.

Under the DO representation, the main SPDE (6)–(8) can be reformulated into the following system of equations [26]:

$$\frac{\partial \overline{u}\left(x,t\right)}{\partial t}=E\left[ \mathcal{L}\left[u\left(x,t;\omega \right);\omega \right]\right].$$

(14)

$$\frac{\partial u_i\left(x,t\right)}{\partial t}={{\displaystyle \prod }}_{V_s^{\perp }}\left[E\left[\mathcal{L}\left[u\left(.,t;\omega \right);\omega \right]Y_j\left(t;\omega \right)\right]\right] C_{Y_i\left(t\right)Y_j\left(t\right).}^{-1}$$

(15)

$$\frac{d{Y}_i\left(t;\omega \right)}{dt}=\langle \mathcal{L}\left[u\left(x,t;\omega \right);\omega \right]-E\left[ \mathcal{L}\left[u\left(x,t;\omega \right);\omega \right]\right],u_i\left(x,t\right)\rangle ,$$

(16)

where ${{\displaystyle \prod }}_{V_s^{\perp }}\left[F\left(x\right)\right]=F\left(x\right)-{{\displaystyle \prod }}_{V_s}\left[F\left(x\right)\right]=F\left(x\right)-F\left(.\right),u_k\left(.,t\right) u_k$ is the orthogonal complement projection, and $C_{Y_i\left(t\right)Y_j\left(t\right)}=E\left[Y_i\left(t;\omega \right)Y_j\left(t;\omega \right)\right]$ is the covariance operator for the stochastic coefficients. The related boundary values are given by:

$$B\left[\overline{u}\left(\xi ,t;\omega \right)\right]=E\left[h\left(\xi ,t;\omega \right)\right].$$

(17)

$$B\left[u_i\left(\xi ,t\right)\right]=E\left[Y_i\left(t;\omega \right)h\left(\xi ,t;\omega \right)\right] C_{Y_i\left(t\right)Y_j\left(t\right).}^{-1}$$

(18)

where $\xi \in \partial D$ and the initial conditions at time $t_0$ have the form:

$$\overline{u}\left(x,t_0\right)={\overline{u}}_0\left(.\right)\equiv E\left[u_0\left(x;\omega \right)\right].$$

(19)

$$u_i\left(x,t_0\right)=u_{i0}\left(x\right).$$

(20)

$$Y_i\left(t_0;\omega \right)=\langle u_0\left(.;\omega \right)-{\overline{u}}_0\left(.\right),u_{i0}\left(.\right)\rangle ,$$

(21)

for $i=1,\dots ,N$, and $u_{i0}\left(x\right)$ are the eigenfunctions associated with the covariance operator $C_{u,u}$.

2.2. Numerical Solution

The DO evolution Equations (14)–(16) include numerical integration for both random and physical spaces. We denote the weights and collocation points of the spatial space by ${\left(x_p,w_p\right)}_{p=1}^{N_s}$ and for the random space by ${\left({\xi }_j,{\gamma }_j\right)}_{j=1}^{N_r}$. We can use Fourier collocation points for $x_p, p=1,\dots ,N_s,$ and Legendre-Gauss collocation for ${\xi }_j, j=1,\dots , N_r$. To discretize time, explicit procedures can be used such as Euler’s method and 4th-order Runge-Kutta (RK4) method. The inner products included in the KL equations can be evaluated as:

• The inner product defined in the spatial domain:

$$\langle p\left(x,t;\omega \right),q\left(x,t;\omega \right)\rangle \approx {\displaystyle \sum }_{p=1}^{N_s}p\left(x_p,t;\omega \right) q\left(x_p,t;\omega \right)w_p.$$

(22)

• The mean value operator defined in random space:

$$E\left[p\left(x,t;\omega \right),q\left(x,t;\omega \right)\right]\approx {\displaystyle \sum }_{j=1}^{N_r}p\left(x,t;{\xi }_j\right) q\left(x,t;{\xi }_j\right){\gamma }_j.$$

(23)

The evolution equations of the DO representation (14)–(16) can be alternatively reformulated into a matrix form as follows:

$$\begin{gathered} \frac{\partial \overline{u}}{\partial t}=E\left[\mathcal{L}\left[u\right]\right]. \\ \frac{\partial u}{\partial t}=\left(p-uG\right)C^{-1}. \\ \frac{dY}{dt}=h. \end{gathered}$$

(24)

where:

$$\begin{gathered} C_{ij}=E\left[Y_i Y_j\right]. \\ {h}_j=\langle \mathcal{L}\left[u\right]-E\left[\mathcal{L}\left[u\right]\right],u_j\left(.,t\right)\rangle . \\ {p}_j=E\left[\mathcal{L}\left[u\right] Y_j\right]. \\ G_{ij}=\langle E\left[\mathcal{L}\left[u\right] Y_j\right],u_i\rangle . \end{gathered}$$

(25)

Practically, in several cases, the initial conditions are deterministic, and the stochasticity emanates from other causes such as random forcing and random coefficients. In this case, the stochastic basis $Y_i$ will be initially zero, and hence the covariance matrix for the stochastic coefficients $C=E\left[Y^{T}Y\right]$ will be singular. For this reason, a hybrid approach, combining the PCM and DO methods, is proposed to overcome this problem [28]. We can start the computation with the PCM method to allow the development of randomness up to a certain converting (switching) time $t_s$, then transform to the DO method and use KL expansion to initialize the DO components $\overline{u} ,\left\{Y_i\right\}$ and $\left\{u_i\right\}$ [33,34].

3. Demonstrating Examples

In this section, we shall consider using the DO and PCM-DO methods to analyze Burgers’ viscous and inviscid models with random parameters. For the numerical computations, we shall consider the following parameters:

$$∆t=0.001,t_f=5, N_s=128, N_r=64, \sigma =1,\overline{v}=1,$$

where $N_s$ refers to the number of collocation points in the physical domain and $N_r$ is the number of points in the random domain. The physical domain is discretized using the Fourier collocation points, while the random space is discretized using the Legendre-Gauss collocation. The RK4 method is employed for the time-integrator. At $t=0$, the stochastic variations are zero; i.e., $Y_i$ and $u_i$ cannot be declared. The DO representation could not be initialized at $t=0$. Therefore, we start the simulation at $t=∆t$ to avoid the singularity issue due to the presence of deterministic (nonrandom) initial conditions. The mean, variance, and stochastic and spatial bases of the DO representation are initiated with the modes obtained from KL expansion.

Alternatively, the hybrid PCM-DO approach is used to overcome the singularity due to deterministic initial conditions. In this case, we start the numerical computation using the PCM method up to a switching time $t_s$ to guarantee that the stochasticity is evolved, then we use the KL decomposition to initialize the DO components and hence switch to the DO method.

For comparison and validation, a reference solution using the PCM method with RK4 is considered. The reference solution will use similar parameters as in the DO and PCM-DO methods.

3.1. Stochastic Burgers’ Equation

Considering the viscous Burgers’ equation with stochastic forcing term in the form [28,35]:

$$\begin{gathered} \frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}=v \frac{{\partial }^{2}u}{\partial x^2}+\frac{1+\xi }{2}\sin \left(2\pi t\right),\hspace{1em}x\in \left[0,2\pi \right] and t\in \left[0, t_f\right], \\ u\left(x,0\right)=g\left(x\right),\hspace{1em}x\in \left[0,2\pi \right], \end{gathered}$$

(26)

where $\xi ~U \left[-\mathit{1},\mathit{1}\right],$ the diffusion coefficient $v=0.04$ with periodic boundary conditions and the initial condition $g\left(x\right)$ taken as:

$$g\left(x\right)=\frac{1}{2}\left(\exp \left(\cos \left(2x\right)\right)-1.5\right) \sin \left(3x+0.74\pi \right).$$

(27)

Apply DO method to obtain the evolution operator $\mathcal{L}$ in the form:

$$\frac{\partial u}{\partial t}=\mathcal{L}\left[u\left(x,t;\omega \right)\right]=-u{u}_x+v{u}_{xx}+\frac{1+\xi }{2}\sin \left(2\pi t\right).$$

(28)

By substituting with KL expansion of $u\left(x,t;\omega \right)$, we get:

$$\mathcal{L}\left[u\left(x,t;\omega \right)\right]=-\overline{u}{\overline{u}}_x-Y_i\frac{\partial }{\partial x}\left(u_i\overline{u}\right)-Y_i{Y}_j{u}_i\frac{\partial u_j}{\partial x}+v\left({\overline{u}}_{xx}+Y_i\frac{{\partial }^2{u}_i}{\partial x^2}\right)+\frac{1+\xi }{2}\sin \left(2\pi t\right).$$

(29)

Apply the mean value operator to get:

$$E\left[\mathcal{L}\left(u\right)\right]=-\overline{u}{\overline{u}}_x-C_{ij}{u}_i\frac{\partial u_j}{\partial x}+v {\overline{u}}_{xx}+0.5\sin \left(2\pi t\right).$$

(30)

where $C_{ij}=E\left[Y_i{Y}_j\right]$ and $Y_i\left(t;\omega \right)$ are zero mean stochastic processes. Multiplying the evolution Equation (26) by $Y_j$, we get

$$\mathcal{L}\left(u\right) Y_j=Y_j\left[-\overline{u}{\overline{u}}_x-Y_i{u}_i\frac{\partial \overline{u}}{\partial x}-Y_k\frac{\partial u_k}{\partial x}\overline{u}-Y_i{Y}_j{u}_i\frac{\partial u_k}{\partial x}\right]+v{Y}_j{\overline{u}}_{xx}+v{Y}_i{Y}_j\frac{{\partial }^2{u}_i}{\partial x^2}+\frac{1}{2}{Y}_j\sin \left(2\pi t\right)+\frac{\xi }{2}{Y}_j\sin \left(2\pi t\right),$$

Apply the expectation operator to get:

$$E\left[\mathcal{L}\left(u\right) Y_j\right]=-\left(C_{ij}{u}_i {\overline{u}}_x+C_{kj }\frac{\partial u_k}{\partial x}\overline{u}+C_{ikj} \frac{\partial u_k}{\partial x} \right)+v C_{ij} \frac{{\partial }^2{u}_i}{\partial x^2}+E\left[ \frac{\xi }{2}{Y}_j\right]\sin \left(2\pi t\right).$$

(31)

where $C_{ikj}=E\left[Y_i{Y}_k{Y}_j\right]$ and $C_{kj}=E\left[Y_k{Y}_j\right]$. After we get $E\left[\mathcal{L}\left(u\right)\right]$ and $E\left[\mathcal{L}\left(u\right) Y_j\right]$, we can apply the matrix form in Equation (25). Then calculate the matrices $h_j, p_j$ and $G_{ij}$. Equation (31) is complicated due to the expectation $E\left[\mathcal{L}\left(u\right) Y_j\right]$ that contains the third moment of the stochastic basis. A deterministic initial condition makes the stochastic coefficient $Y_i, i=1,\dots N$ be zero initially, and hence the covariance matrix $C$ for $Y_i$ becomes singular. The hybrid PCM-DO method will be used to avoid this issue. For the numerical solution, we shall use the parameters $N_r$, $N_s$, $∆t$, and discretization methods as shown above.

In Figure 1, the mean and variance are shown at $t_f=5$ and $N=6$. They agree well with the reference solution. The $L_2$ error for the mean and variance for different values of modes are shown in Figure 2. We present errors in mean and variance for three values of modes: $N=3$, $N=5$, and $N=7$. The solution at $t_s=2$ and $N=7$ is the best among the three cases. To this end, it can be noticed that the accuracy of the solution is affected by lower modes.

Figure 3 shows the high convergence rate detected with the number of modes. The $L_2$ error reduces as the number of modes increases. This demonstrates that the DO decomposition has a high convergence rate for this nonlinear case.

Figure 4 shows that the hybrid approach is affected by the switching time and number of modes. We observe that $L_2$ error for the mean and variance at $t_s=1$ and $N=6$ is the best among other values. For $t_s=2$ and $N=6$, the solution has the same accuracy but takes more computational time. The computational time is compared when using the DO and PCM for this test case using different spatial nodes N_r. The DO is more efficient, with less computational time than the PCM, as shown in Figure 5.

Figure 6 shows the eigenvalues variation as the solution evolves till $t_f=5$. The computed eigenvalues are compared with the reference solution. A good agreement is noticed as shown in the figure. It is observed that the lower eigenvalues ${\lambda }_8 \mathrm{and} {\lambda }_9$ from the mode $N=9$ have an order: ${10}^{-15} \mathrm{and} {10}^{-16}$, respectively, making the covariance matrix extremely ill-conditioned; it becomes singular and hence sensitive to any perturbations.

3.2. Quasilinear Inviscid Burgers’ Equation

Consider the stochastic quasilinear inviscid Burgers’ equation with a random forcing in the form [36]:

$$\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}=\frac{1+\xi }{2}\sin \left(2\pi t\right);\hspace{1em}u\left(x,0\right)=g\left(x\right),$$

(32)

where $x\in \left[0,2\pi \right] , t\in \left[0, t_f\right] ,$ $\xi ~ U\left[-1,1\right]$ with periodic boundary conditions and the initial condition is:

$$g\left(x\right)=\frac{1}{7} \sin x.$$

(33)

Applying the DO decomposition to get the evolution operator $\mathcal{L}$ along with other forms:

$$\mathcal{L}\left[u\left(x,t;\omega \right)\right]=\frac{\partial u}{\partial t}=-u{u}_x+\frac{1+\xi }{2}\sin \left(2\pi t\right).$$

(34)

$$E\left[\mathcal{L}\left(u\right)\right]=-\overline{u}{\overline{u}}_x-C_{ij}{u}_i\frac{\partial u_j}{\partial x}+0.5\sin \left(2\pi t\right).$$

(35)

$$E\left[\mathcal{L}\left(u\right) Y_j\right]=-\left(C_{ij}{u}_i {\overline{u}}_x+C_{kj }\frac{\partial u_k}{\partial x}\overline{u}+C_{ikj} \frac{\partial u_k}{\partial x} \right)+E\left[ \frac{\xi }{2}{Y}_j\right]\sin \left(2\pi t\right),$$

(36)

where $C_{ikj}=E\left[Y_i{Y}_k{Y}_j\right]$. We shall use similar parameters as in the above example in Section 3.1.

Figure 7 demonstrates that errors in the mean and its variance are identical for both the DO and the PCM methods. We present mean and variance at time $t_f=6$ as shown in Figure 8. They have a good agreement compared with the reference solution. Figure 9 shows the performance of the mean solution $\overline{u}\left(x,t\right)$ using the DO method as the number of time iterations increases along with the contour lines. The solution is illustrated at $t_f=6$ at switching time $t_s=1$ and number of modes $N=4$ with time increment $∆t={10}^{-3}$. We observe that the wave shape tends to be sharper with time. The $L_2$ error of the DO method for different modes number $N$ are shown in Figure 10. We present the errors in the variance and mean for three modes $N=2,N=3 \mathrm{and} N=4$ with the same switching time $t_s=1$. As shown, and as expected, the error due to mode $N=4$ is less than errors due to other modes. Therefore, as the number of modes increases, the error is reduced. Figure 11 shows the high rate of error convergence with increasing the modes number. This implies that the DO representation is also an efficient choice for this nonlinear SPDE.

As it is shown in Figure 12, the $L_2$ error of the solution is computed for three switching times $t_s=0.5, t_s=1$ and $t_s=2$ with the same number of modes $N=4$. We can observe that the error for the mean at $t_s=2$ is better than the errors due to other switching times. After that, the error becomes the same for the three modes. On the other hand, the error in variance is approximately the same for all switching times, so we can use $t_s=1$ as it has a smaller computational time.

Figure 13 shows the time growth of eigenvalues till $t_f=6$. The computed eigenvalues are compared with those from the reference solution and a good agreement is achieved. It is observed that the smaller eigenvalues ${\lambda }_3$ and ${\lambda }_4$ have an order ${10}^{-15} \mathrm{and} {10}^{-12}$, respectively. This means the covariance matrix is highly ill-conditioned and becomes singular as the modes number increases.

4. Shock Wave Occurrence

Consider the stochastic inviscid Burger’s equation with random forcing as follows:

$$\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}=\frac{1+\xi }{2}\sin \left(2\pi t\right);\hspace{1em}u\left(x,0\right)=g\left(x\right),$$

(37)

where $x\in \left[0,2\pi \right] , t\in \left[0, t_f\right] ,$ $\xi ~ U\left[-1,1\right]$ with periodic boundary conditions and the initial condition to be taken as:

$$g\left(x\right)=\frac{1}{3} sinx.$$

(38)

This model is the same model given in (32), but with different initial condition $g\left(x\right)$ that allows for shock formation. We consider the solution up to a final time $t_f=4$. The above described DO method in Section 2 is applied to the model Equation (37). The convective term $u \partial u/\partial x$ should be handled using any suitable technique such as first and/or second-order upwind scheme. This will prevent oscillatory behavior of the solution, especially near the shock wave. Figure 14 shows the solution mean and variance at the occurrence of the shock wave. The solution is illustrated at switching time $t_s=1$, and the number of modes considered is $N=4$. As the waveform steepens, the edges form causing a shock wave. We can notice how variance become relatively large close to the shock wave. The maximum value of variance occurs exactly at the shock wave. The mean field $\overline{u}\left(x,t\right)$ of the solution is shown in Figure 15. We present the mean field and the contour lines’ evolution up to the breaking time. The sinusoidal wave shape sharpens as the time evolves till the shock wave occurs. The mean and variance for the solution is sketched at different values of time as shown in Figure 16. Time values $t=1.1, t=2.5, t=3.5,$ and $t=4$ are considered in Figure 16. It can be observed that the sinusoidal wave shape is distorted. The distortion of the wave profile starts from the initial data, and the dispersion of the wave profile evolves with the time increase.

5. Conclusions

The Burgers’ model with stochastic terms is an important problem which appears in many applications such as the flow turbulence. Due to lack of analytical methods, numerical techniques are required to analyze such important models. The stochastic Burgers’ model is studied using DO and PCM methods in the current study. The DO method introduces a computational low-dimensional framework that enables to capture the relevant system dynamics for the stochastic Burgers’ model. The PCM is combined with the DO to avoid the singularity of the DO evolution in case of deterministic initial conditions. The approach is based on solving the model using the PCM up to a certain switching time at which the stochasticity occurs, and then switching to the DO method. The two strategies, DO and PCM-DO, are successfully applied to the stochastic Burgers’ model and are compared in terms of performance and accuracy with a reference solution. The stochastic viscous and inviscid Burgers’ models combined with random forces are considered with the two methods. We derived the quantities of the stochastic coefficients and spatial bases for these models. When analyzing the test models, the DO method shows a high level of agreement with the reference solution and has a convergence rate increase with the number of modes.

Overall, the DO method is as accurate as the PCM method, but with higher efficiency. This explains why the DO method could be a viable and gainful replacement for the PCM. We also analyzed the shock wave formation using the DO by considering the inviscid Burgers’ equation with suitable discretization of the convective tern to avoid oscillatory solution. The sinusoidal wave shape sharpens and distorts as the time evolves till the shock wave occurs. The mean solution shows the shock occurrence while the variance explains the magnification of oscillations close to shock formation. Further research should focus on more systematic studies to determine the suitable switching time from PCM to DO methods. Moreover, other applications are to be considered using the current presented techniques.