Abstract

There are some previous works on designing efficient and high-order numerical methods of density estimation for stochastic partial differential equation (SPDE) driven by multivariate Gaussian random variables. They mostly focus on proposing numerical methods of density estimation for SPDE with independent random variables and rarely research density estimation for SPDE is driven by multivariate Gaussian random variables. In this paper, we propose a high-order algorithm of gPC-based density estimation where SPDE driven by multivariate Gaussian random variables. Our main techniques are (1) we build a new multivariate orthogonal basis by adopting the Gauss–Schmidt orthogonalization; (2) with the newly constructed orthogonal basis in hand, we first assume the unknown function in the SPDE has the stochastic general polynomial chaos (gPC) expansion, second implement the stochastic gPC expansion for the SPDE in the multivariate Gaussian measure space, and third we obtain and numerical calculation deterministic differential equations for the coefficients of the expansion; (3) we used high-order algorithm of gPC-based for density estimation and moment estimation. We apply the newly proposed numerical method to a known random function, stochastic 1D wave equation, and stochastic 2D Schnakenberg model, respectively. All the presented stochastic equations are driven by bivariate Gaussian random variables. The efficiency is compared with the Monte-Carlo method based on the known random function.

1. Introduction

In modern science, no matter numerical model or physical model, the study of uncertainty (randomness) is more important and concerned. Such uncertainties arise from one or more parameters of the model, sometimes related to each other. Because of such uncertain parameters, the solution of the model is also uncertain. For example, in the field of biochemistry [1–3], in the field of structural engineering [3, 4], and in hydrology engineering [3, 5].

For the random solution, there are usually two main research aspects of moment estimation and density estimation from the perspective of statistical properties. The research of moment estimation, such as mean (expectation) and variance, is based on obtaining the numerical solution of the model by numerical calculation. The most direct method to obtain numerical solutions is the Monte-Carlo method, which takes samples from the sample space of parameters describing randomness, assigns values to random parameters, and numerically calculates the random solutions. However, this method has a large amount of calculation and a low convergence rate [6]. If high precision approximation is required, stochastic general polynomial chaos (gPC) expansion is a preferred choice, such as the stochastic Galerkin spectral method. We also can consider the stochastic perturbation-based finite element method, and this method is a very efficient alternative for higher-order spectral methods. The theoretical basis of the method is Taylor’s expansion of all uncertain parameters and state functions [7].

For the stochastic Galerkin spectral method, how to use the weight function of random variables (joint probability density function) to find the orthogonal basis is a major challenge. Especially when dealing with multivariate random variables, a multivariate orthogonal basis is difficult to find [8].

Orthogonal basis can be used in relevant methods, including measure-consistent chaos expansion [9, 10], domination method [11], Gauss–Schmidt orthogonalization [12–18], and mapping transformation method [8, 11, 19, 20]. After obtaining the approximation of the random function, the moment estimation can be carried out, and the mean or variance of the random solution can be estimated as a whole.

In this paper, we discuss the differential equations driven by multivariate Gaussian random variables. On the basis of moment estimation, we estimate the density of random solution. When the random variable is univariate, standard kernel density estimation [3, 21, 22] can be used from a statistical point of view. In the study of Ditkowski et al. [3], numerical solutions of random functions and density of stochastic differential equations driven by univariate random variables and multivariate linear independent random variables are discussed, mainly gPC-based estimation and Spline-based estimation. For the multivariate random variable (DEPENDENT) situation, it introduced briefly the application of multivariate-spline-based estimation [3].

Based on the above excellent results and methods, this paper aims at differential equations driven by multivariate random Gaussian variables. A high-order algorithm of gPC-based for density estimation based on polynomial chaos expansion is designed to estimate the moment and density of the random solution [3]. Our main techniques are (1) we build a new multivariate orthogonal basis by adopting the Gauss–Schmidt orthogonalization; (2) with the newly constructed orthogonal basis in hand, we first assume the unknown function in the stochastic partial differential equation (SPDE) has the stochastic gPC expansion, second implement the stochastic gPC expansion for the SPDE in the multivariate Gaussian measure space, and third we obtain and numerical calculation deterministic differential equations for the coefficients of the expansion; (3) we used high-order algorithm of gPC-based for density estimation and moment estimation. We apply the newly proposed numerical method to a known random function, stochastic 1D wave equation, and stochastic 2D Schnakenberg model, respectively. The efficiency is compared with the Monte-Carlo method based on the known random function.

This article is organized as follows: in Section 2, we discuss the high-order stochastic chaos expansion of random functions and introduce the stochastic Galerkin spectral method and the high-order approximation of the expansion coefficients. In Section 3, the high-order algorithm of gPC-based for density estimation based on polynomial chaos expansion is introduced, and the numerical calculation process of moment estimation and density estimation is introduced. In Section 4, the moment and density estimation of a known random function are calculated numerically, and the efficiency problem is compared with the Monte-Carlo method. The proposed algorithm is applied to stochastic 1D wave equation and stochastic 2D Schnakenberg model. In Section 5, the algorithm and numerical results proposed in this paper are summarized.

2. High-Order Stochastic gPC Expansion

We consider the boundary value problem of SPDE driven by multivariate Gaussian measurewhere the unknown function . is the spatial variable, are some differential operator with respect to . is a K-dimensional random vector with Gaussian measure. are some known function. is the -dimensional vector of random input variables, whose Joint probability density function (). typically represents the uncertainties in the model [8, 23].

2.1. Stochastic General Polynomial Chaos Expansion

Given a finite nonnegative integer and , we define a truncated multi-indices set as follows [24]:

The polynomial chaos expansion for random input functions , , and and the unknown function [8]where the polynomial chaos orthogonal basis can be constructed with Gauss–Schmidt orthogonalization method (more details, see [18]). Then, we can obtain the following finite-dimensional and deterministic PDE problem: for all , finding such that

2.2. Higher-Order Calculation of Coefficients

In order to high-order to solve Equation (4), the calculation of coefficients as follows:Denote the integral transformation for where . Then Equation (5) can be rewritten as follows:

We choose mesh sizes with for even positive integers, and let the grid points be as follows:

Then Equation (7) can be rewritten as follows where used trapezoidal rules, denotes

In order to high-order approximation, the length of every dimension needs to a lot more than usual.

3. Density and Moment Estimation

For random function , we may wish to know its statistical properties as follows:

1. Moment estimation.

2. Density estimation.where is the cumulative distribution of , .

3.1. Higher-Order Algorithm of gPC-Based

3.2. Accuracy of Algorithm

The error of density estimation comes from the approximation for . Denote and the set equal to .

For some constant [3, 25–27].where represent standard deviation. If is analytic, the truncated expansion (Equation (3)) has the exponential accuracy [3, 27, 28].

4. Numerical Experiments

In this section, first, we text the accuracy of Algorithm 1 in moment estimation and density estimation for a known random function driven by bivariate Gaussian random variables. Secondly, let’s apply Algorithm 1 for the stochastic 1D wave equation and stochastic 2D Schnakenberg model with the same random variables as the known random function.

The bivariate Gaussian random variables (i.e., they are bivariate Gaussian random variables). are some known constants ( is taken in our calculation later).

The covariance of is .

The joint probability density of the bivariate Gaussian distribution is as follows:

4.1. Example

4.1.1. gPC Chaos

We approximate a known random function with the truncated polynomial chaos expansion. Let the random function be defined as follows:

Let , then . We numerically computed multivariate orthogonal polynomial basis based on Gauss–Schmidt orthogonalization method [18], the shows as Table 1.

4.1.2. Moment Estimation

Furthermore, we use the SG method to find the coefficients of random function (Equation (19)), and Table 2 summarizes our result.

Table 3 shows the exact mean and the approximate mean of random function (Equation (19)), the exact variance and the approximate variance, and the error between them, respectively.

4.1.3. Density Estimation

Based on Algorithm 1 and Tables 1 and 2, we can obtain the truncated expansion for the random function and calculated its probability density function (PDF). In order to sampling, we taken each dimension has the same sample size, .

are sample points according to . Here, approximation uses sample points. We used square Euclidean distance to describe the error of gPC-based density estimation.

Figure 1(c) shows the error of gPC-based are larger (red dotted line) when the number of basis function , and the error has the spectral accuracy when (blue dotted line).

(a)

(b)

(c)

(a)
(b)
(c)

Figure 1 

(a) Contours of the random function . (b) The PDF of (black diamond dotted line), the red line and blue line represent its approximation by multidimensional spline-based and gPC-based Algorithm 1, respectively. (c) The error of PDF by gPC-based approximation with the number of basis function.

4.1.4. Efficiency

In this section, we focus on moment estimation and density estimation for random function in contrast to Monte-Carlo simulation with different sample sizes.

Table 4 shows the mean and variance of by Monte-Carlo simulation; the error of mean and variance gets smaller and smaller when the sample size of simulation multiplies. In contrast to Table 3, it can be seen that high-order methods are spectral accuracy when the degrees of basis function are two orders (). The error of mean and variance by Monte-Carlo simulation can be spectral accuracy when the sample size is enormous amount ().

Table 5 shows the PDF error of by Monte-Carlo simulation. In contrast to Figure 1, the proposed method can achieve high precision with only six expansion terms, and Monte-Carlo simulation requires a huge amount of computation to achieve close precision. Although the accuracy of Monte-Carlo simulation is limited and the calculation cost is high, this method is easy to implement and is very simple and useful. Then it can be seen that high-order methods are superior to the Monte-Carlo simulation, as it gives sufficiently accurate results [29, 30].

4.2. Application

In this section, we used the gPC-based density estimation to the stochastic 1D wave equation and Stochastic 2D Schnakenberg model. We assume that the unknown function has periodic boundary conditions in the or direction, and we will use the Fourier spectral method to approximate the unknown function in the or direction [31–33]. The are bivariate Gaussian distribution, same as Equation (19).

4.2.1. Stochastic 1D Wave Equation

The stochastic 1D wave equation.where .

Table 6 shows the mean and variance of at different time when . With the passage of time, the mean amplitude of the wave gradually decreases, and the variance gradually decreases and tends to zero, while Figure 2 (top and middle) also demonstrates this phenomenon. When , the evaluated at point 0.5, so, the value of PDF for are huge at location u nearby point 0.5.

Figure 2 

The mean amplitude of Equation (23) solution at different times (t = 0, 5, 10, 20 (top)). The trend fluctuations of mean amplitude for Equation (23) solution at every time (middle). The PDF of Equation (23) solution at different times (t = 0, 5, 10, 20 (bottom)).

As the wave energy dissipates, the value of PDF for and gradually close to the location of u nearby zero.

Until the wave energy almost runs out, the value of PDF for evaluated at location of u nearby zero and the interval length of is small than before. At the location of , the value of PDF for almost equal to 1.where is a small neighborhood with centered on 0 point.

4.2.2. Stochastic 2D Schnakenberg Model

The Stochastic 2D Schnakenberg model and Schnakenberg model belong to the reaction-diffusion system.where , .

In stochastic 2D Schnakenberg model, and can be described as the concentration of some material. Table 7 shows the mean and variance of at different time when . From Table 5, we know, as the beginning, the mean and variance of is zero, and the mean and variance of is increase when time increase. when , the mean and variance of almost same as and , just a small change when time increase.

When , the equal zero, so, the value of PDF for evaluated at location of u nearby zero and the interval length of is small.where is a small neighborhood with centered on 0 point.