## Abstract

Approximating nonlinearities in stochastic partial differential equations (SPDEs) via the Wick product has often been used in quantum field theory and stochastic analysis. The main benefit is simplification of the equations but at the expense of introducing modelling errors. In this paper, we study the accuracy and computational efficiency of Wick-type approximations to SPDEs and demonstrate that the Wick propagator, i.e. the system of equations for the coefficients of the polynomial chaos expansion of the solution, has a sparse lower triangular structure that is seemingly universal, i.e. independent of the type of noise. We also introduce new higher-order stochastic approximations via Wick–Malliavin series expansions for Gaussian and uniformly distributed noises, and demonstrate convergence as the number of expansion terms increases. Our results are for diffusion, Burgers and Navier–Stokes equations, but the same approach can be readily adopted for other nonlinear SPDEs and more general noises.

## 1. Introduction

The Wick product was originally introduced by Wick [1] in quantum field theory to reduce a product of creation and annihilation operators to sums of products of pairs of these operators. A similar concept was introduced in stochastic analysis by Hida & Ikeda [2] to analyse Hilbert spaces with reproducing kernels arising from multiple Wiener integrals. Since then, the Wick product has played an increasingly important role in many fundamental problems involving stochastic partial differential equations (SPDEs) [3,4] and in the theory of stochastic integration, as it can be considered as a generalization of the Skorokhod–Malliavin integral [5]. Replacement of product-type nonlinearities in SPDEs by Wick products can also be viewed as a stochastic approximation technique formulated at the level of the equation, i.e. before any discretization. For example, consider the stochastic Burgers equation , driven by Gaussian space–time white noise . The Wick formulation of this problem is [6], where denotes the Wick product operator.^{1} Similarly, one can define the Wick approximation of other nonlinear SPDEs, for example, the diffusion equation with random inhomogeneous diffusivity recently studied by Wan *et al.* [8].

In this paper, we present theoretical and numerical results for a new type of approximation, i.e. the Wick–Malliavin approximation. The key idea relies on expanding the nonlinearities appearing in SPDEs in a Taylor-like series involving Wick products and Malliavin derivatives [9,10]. In this generalized framework, Wick-type SPDEs can be regarded as a zero-order truncation of the Wick–Malliavin series. Higher-order approximations can be constructed in a systematic way by including terms involving Malliavin derivatives of higher order.

The Wick approximation to SPDEs has interesting properties. In particular, it represents an *unbiased perturbation* of the SPDE, which leads to significant computational advantages for either white or coloured noise. In fact, as we will see, the Wick propagator, i.e. the system of equations for the deterministic coefficients of the polynomial chaos expansion of the solution, is *quasi-linear* and *sparse* for Wick-type equations, including equations with polynomial Wick nonlinearities [6,9], Wick-multiplicative Gaussian noise [4,8,11], as well as noises represented by nonlinear transformations of Gaussian processes. By ‘quasi-linear’ here, we mean that only the equation for the mean field is nonlinear, whereas the rest of the system is linear, forming a hierarchy that can be readily solved. In particular, the propagators of the Wick approximation to the stochastic Burgers or Navier–Stokes (NS) equations are *lower triangular* systems. By contrast, the propagators of the standard stochastic Burgers or NS equations are infinite-dimensional systems of fully coupled nonlinear PDEs.

On the other hand, the Wick approximation is a rather draconian stochastic truncation of the nonlinearities appearing in SPDEs, and hence we need to address the question of accuracy, that is how far the Wick stochastic solution deviates from the stochastic solution of the original SPDEs. This issue has been answered only partially in the past, for example, by Wan *et al.* [8], where it was shown numerically that for elliptic SPDEs, the variance of the Wick solution was close to the stochastic solution for small perturbations but deviated greatly for large noise levels. Here, we revisit this issue and also study higher-order Wick–Malliavin approximations to the stochastic Burgers and NS equations. As we will see, the corresponding propagators are no longer lower triangular as in the Wick setting. Specifically, the coupling strength between the chaos modes, i.e. the number of terms above the principal diagonal of the propagator, is related to the truncation order in the Wick–Malliavin series. A low-order truncation yields a quasi-lower triangular system that can be efficiently solved by using operator-splitting techniques. A fundamental question is whether the Wick–Malliavin approach is limited to Gaussian perturbations or if it can be extended to other types of noises. It turns out that with the appropriate definition of Wick- and Malliavin-type derivatives, both the theory, and the numerical methodology can be generalized to *non-Gaussian* processes. In particular, we have studied Gaussian and uniformly distributed perturbations of Burgers and NS equations, but the same approach can be readily extended to other SPDEs with more general noises.

This paper is organized as follows. In §2, we introduce the Wick–Malliavin series expansion for random fields, and discuss the coupling structure and the sparsity of the propagator associated with quadratic nonlinearities. In §§3–5, we present the application of the Wick–Malliavin series to the diffusion equation with inhomogeneous random diffusivity, the Burgers equations driven by random noise and the NS equations subject to random boundary and random initial conditions, respectively. The accuracy and the computational efficiency of the Wick–Malliavin approximation is discussed in §6. Finally, the main findings and their implications are summarized in §7. We also include a brief appendix dealing with Wick–Malliavin rules for univariate Hermite and Legendre polynomials.

## 2. Wick–Malliavin series expansion

It has been recently shown by Mikulevicius & Rozovskii [9] that a random polynomial nonlinearity, such as the nonlinear advection term in the NS equations, can be expanded in a Taylor-like series involving *Wick products* [7] and *Malliavin derivatives* [10]. Specifically, given two square-integrable random fields *u* and *v*, it has been shown that
2.1where denotes the Wick product and is the Malliavin derivative of order *p*. The series (2.1) converges with probability one for sufficiently smooth random fields. It can also be used as a starting point for constructing approximations of *finite stochastic order*. For example, the nonlinear advection term in the NS equations can be approximated as
2.2where *Q* denotes the highest stochastic order in the Wick–Malliavin expansion. In particular, the zero-order approximation is known as the *Wick approximation*. Similarly, the first-order (or level-one) Wick–Malliavin approximation is .

### (a) Wiener–Hermite representations

Let us consider a square-integrable scalar random field *u* depending on *M* (finite or infinite) independent identically distributed (i.i.d.) Gaussian random variables (*ξ*_{1},…,*ξ*_{M}). Denote by the set multi-indices ** k**=(

*k*

_{1},…,

*k*

_{M}) such that and 2.3By the Cameron & Martin [12] theorem (see also [13]),

*u*admits the following

*L*

_{2}-convergent expansion 2.4are multivariate monic Hermite polynomials of Gaussian random variables (

*h*

_{kj}(

*ξ*

_{j}) are univariate monic Hermite polynomials of order

*k*

_{j}). The quantities are deterministic fields (chaos modes). The Wick product between two elements

*h*

_{k}and

*h*

_{i}can be defined as (see [7,14]) 2.5The Malliavin derivative of

*h*

_{k}can be defined as 2.6where

*q*_{j}is a multi-index such that only the

*j*th component is equal to one, whereas all the others are zero. The definition we give here is for monic Hermite polynomials of Gaussian random variables (see [11]). Mikulevicius & Rozovskii [9] and Nualart [10] give a much more general definition for orthogonal polynomial functionals of generalized processes, such as cylindrical Brownian motion. In the latter case, the Malliavin derivative involves a representation of the random space in terms of a complete orthonormal basis. The

*p*th order derivative, , can be computed by induction from (2.6). An important result, recently proved by Mikulevicius & Rozovskii [9], is the following.

### Proposition 2.1

*Let h*_{i} *and h*_{j} *be elements of the Cameron–Martin basis. Then, with probability one*,
2.7

By linearity, the statement of the proposition allows us to conclude that if *u* and *v* are two square-integrable scalar random fields with Wiener chaos expansions converging sufficiently fast, then equation (2.1) holds true.

The normalized projection of the *Q*th-order truncation of (2.7) on the basis {*h*_{k}} gives us the very important quantity
2.8where 〈·〉 denotes the statistical expectation operator. As we will see shortly, the Wiener chaos projection of the Wick–Malliavin approximation to SPDEs with quadratic nonlinearities—such as the NS equations—yields a propagator whose coupling structure is completely defined by (2.8). Specifically, for each fixed value of the multi-index ** k**, the matrix
2.9gives us the coupling structure of the chaos modes in the

**th equation of the propagator, which ultimately affects the computational cost, efficiency and accuracy of the Wick–Malliavin approximation.**

*k*The structure of the coupling matrices (2.9) is shown in figure 1. Note that the Wick approximation (*Q*=0) is always *lower triangular* and *very sparse* for *M*>2 (see also table 1). On the other hand, the full Wiener chaos system is dense.

### (b) Non-Gaussian noise

By using the linearization theory of orthogonal polynomials [15] and the completeness result proved by Ernst *et al.* [13], one can generalize (2.4)–(2.8) to the non-Gaussian case. In principle, this allows us to define Wick- and Wick–Malliavin-type approximations to nonlinear SPDEs driven by non-Gaussian noise (see §5). In particular, in the appendix, we provide a possible generalization of (2.4)–(2.7) to univariate Legendre polynomials of uniform random variables. The corresponding analytical form of (2.8) is summarized in table 2.

In the following, we present the application of the Wick–Malliavin approximation technique to several fundamental SPDEs of mathematical physics. In particular, we consider stochastic diffusion, Burgers and NS equations.

## 3. Diffusion equation

Let us consider the following boundary-value problem:
3.1where *a*(** x**;

*ω*) is a

*lognormal*homogeneous random field, and the forcing term

*f*is assumed to be deterministic for the sake of simplicity. By using formula (2.1), we expand the direct product between the random fields

*a*(

**;**

*x**ω*) and ∇

*u*(

**;**

*x**ω*) as 3.2where we have truncated the infinite series to the order

*Q*. A substitution of equation (3.2) into equations (3.1) yields the following

*Q*th-

*order Wick–Malliavin approximation*: 3.3In particular, the zero-order truncation (

*Q*=0) of equation (3.3) gives us the Wick approximation 3.4This model has been studied in the past, for example, by Wan

*et al.*[8], Theting [16] and Lototsky

*et al.*[17]. The main motivation was that the Wick product is consistent with the Skorokhod stochastic integral, and it is expected that the solution can be smoothed to some extent by the convolution in the probability space.

### (a) Wick–Malliavin propagator

We consider the following form of the diffusion coefficient:
3.5is a random field expanded in a Karhunen–Loève series, and {*ξ*_{i}} is a set of i.i.d. normal random variables.^{2} The Wiener chaos expansion of (3.5) can be easily obtained as
3.6where |** k**| is defined in (2.3) and . We also consider the Wiener chaos expansion of the solution to equation (3.3),
3.7A substitution of (3.6) and (3.7) into (3.3) and subsequent Galerkin projection on the basis yields the

*Q*th-

*order Wick–Malliavin propagator*3.8where is defined in (2.8). For

*Q*=0 (Wick approximation), the propagator (3.8) becomes

*lower triangular*(figure 1). This means that we can effectively solve each equation of the propagator by using forward substitution. The accuracy of the Wick approximation to stochastic elliptic problems with lognormal coefficients was first studied by Wan

*et al.*[8], where a second-order convergence rate with respect to the noise amplitude

*σ*in (3.5) was established. In this paper, this result is further generalized to Wick–Malliavin approximations. In particular, in §6, we show numerically that each Malliavin derivative added in the propagator improves the convergence order with respect to the noise amplitude

*σ*by an order of two. Specifically, we have convergence order 2 for

*Q*=0 (Wick approximation), 4 for

*Q*=1 (level-one Wick–Malliavin approximation) and 6 for

*Q*=2 (level-two Wick–Malliavin approximation). Remarkably, the same result holds for the randomly forced Burgers equation and, most likely, for more general SPDEs with quadratic nonlinearities. The improvement in the convergence rate comes, however, at the expense of a weak coupling between the equations of the propagator, which is no longer lower triangular as in the Wick setting.

## 4. Burgers equation

Let us consider the following initial/boundary-value problem:
4.1where *ξ*_{k}(*ω*) are i.i.d. normal random variables, , , and *σ*≥0 modulate the amplitude of the forcing term. We expand the nonlinear advection term *u*∂*u*/∂*x* in a truncated Wick–Malliavin series in the form (2.2). This yields the following *Q*th-*order Wick–Malliavin approximation*:
4.2The case *Q*=0 corresponds to the Wick approximation and has been theoretically studied by Kaligotla & Lototsky [6], see also Grothaus *et al.* [3] and Holden *et al.* [4].

### (a) Wick–Malliavin propagator

We represent the solution to equation (4.2) (for arbitrary *Q*) in a Wiener–Hermite series as
4.3where *h*_{k} are multivariate Hermite polynomials and ** k** is a multi-index in a suitable index set . A substitution of equation (4.3) into equation (4.2) and subsequent Galerkin projection onto the space spanned by

*h*

_{k}yields the

*Q*th-

*order Wick–Malliavin propagator*4.4where the coefficients are defined in equation (2.8) (see also figure 1). If we set

*Q*=0 in equations (4.4) and (2.8), and rearrange the double summation, we obtain the

*Wick propagator*4.5This is a lower triangular system that can be readily solved in a sequence by forward substitution; moreover, only the first equation of the system, i.e. the one for

**=**

*k***0**is nonlinear, whereas all the others are linear.

### (b) Stochastic perturbation method

There exists an interesting connection between the Wick approximation and the classical first-order perturbation method in terms of *σ*. In order to exploit such a connection, let us assume that the solution to (4.1) is analytic in the parameter *σ*. This allows us to write the power series expansion as
4.6where *U*_{j}(*x*,*t*;*ω*), *j*=1,2,…, are zero-mean random fields. If we substitute this representation into equations (4.1), equate the terms having the same power of *σ* and truncate the series expansion (4.6) to first-order in *σ*, then we obtain the *linear approximation*
4.7Note that the amplitude of the perturbation *σ* no longer appears in the above system, which means that we can solve it once, obtain the modes *U*_{0} and *U*_{1}, and then construct a solution in the form (4.6) for the desired value of *σ* (provided *σ* is reasonably small). Next, we represent the random field *U*_{1}(*x*,*t*;*ω*) in a (zero-mean) truncated Wiener–Hermite series as
4.8A substitution of equation (4.8) into the system (4.7) and subsequent Galerkin projection onto the space spanned by yields
4.9It is interesting to compare the Wick propagator (4.5) with the system of equations (4.9). The first equation in both systems is the same and it is the only nonlinear equation of the whole set. The structure of the equations for |** k**|=1 is similar apart from a constant factor (

*σ*) that can be rescaled. The remaining part of the system is apparently different. This suggests that for small

*σ*the error between the Wick–Burgers and the full Burgers equation should scale as

*σ*

^{2}because this is the first equation that differs in the systems (4.5) and (4.9). This prediction is confirmed by numerical results in §6. Higher-order perturbation expansions in the form (4.6) can also be considered, leading to systems with lower triangular structure. From a computational viewpoint, however, it is not convenient to go beyond first order because the random fields {

*U*

_{1},

*U*

_{2},…} appearing in (4.6) have to be represented in a Wiener chaos expansion. This means, in particular, that a second-order perturbation expansion doubles the number of unknown chaos modes with respect to a first-order expansion. On the other hand, the level-one Wick–Malliavin approximation keeps the number of unknowns constant, but it introduces a weak coupling between the equations of the propagator, which is no longer lower triangular as in the Wick setting.

## 5. Navier–Stokes equations

Let us consider the incompressible NS equations
5.1subject to a random initial condition or random boundary conditions [18,19]. The nonlinear advection term in (5.1) can be represented via the Wick–Malliavin expansion. To this end, we simply substitute the series (2.2) into equations (5.1), to obtain the following *Wick–Malliavin–Navier–Stokes* (WMNS) equations:
5.2In particular, the zero-order approximation, i.e. the *Wick approximation*, is obtained by setting *Q*=0 in equations (5.2). This yields the *Wick–Navier–Stokes* (WNS) equations
5.3The theoretical properties of equations (5.3) have been recently studied by Mikulevicius & Rozovskii [9]. Let us briefly review some of the main results. First of all, the WNS equations are an *unbiased perturbation* of the deterministic NS equations. In other words, the mean of the solution to equations (5.3) is a solution to the deterministic equations (5.1) with averaged initial and boundary conditions. Another property is that the solution to equations (5.3) exists and is unique but not necessarily square integrable.^{3} As a consequence, the variance of the solution diverges to infinity as time increases. To overcome this problem, Mikulevicius & Rozovskii [9] proposed *renormalization* that uses scaling based on Kondratiev's spaces and Catalan numbers. The reason for the latter is that Catalan numbers arise in the asymptotics of convolutions, and the nonlinear term in equation (5.6) is a convolution. As we will see in the numerical results of §6*c*, the *a priori* rescaling of the WNS solution based on Catalan numbers yields results that are surprisingly very close to the optimal (in the second-order moment sense) rescaling computed on the basis of the solution to the full NS equations.

### (a) Wick–Malliavin propagator

Let us expand the velocity and pressure fields in a finite-dimensional Wiener–Hermite series,
5.4A substitution of equations (5.4) into equations (5.2) and subsequent Galerkin projection onto the space spanned by yields the *Wick–Malliavin propagator*
5.5where the coefficients are defined in equation (2.8) (see also table 2). Setting *Q*=0 in equation (5.5) yields the following *Wick propagator*:
5.6This system is *lower triangular*, *very sparse* for multiple random variables, and it involves only one nonlinear equation, i.e. the one for ** k**=

**0**; the rest of the system is linear. This attractive feature, which can yield significant computational advantages, immediately leads us to the question of accuracy of Wick-type or Wick–Malliavin-type NS equations, relative to the full NS equations. This question is addressed in detail in §6.

It is interesting to note that the propagator of the WNS equations (5.3) has the same structure as the system for the coefficients of *formal power series solutions* to the NS equations (5.1) in the form
5.7where ** ξ**=(

*ξ*

_{1},…,

*ξ*

_{M}). In fact, if we substitute equations (5.7) into equations (5.1) and collect the terms having the same powers of

**, we obtain exactly the system (5.6). The main technical reason for this is that the basis functions**

*ξ**h*

_{k}and

*ξ*^{k}appearing in equations (5.4) and (5.7) have exactly the same type of product structure relative to the Wick and the dot product, respectively. Indeed, by definition, we have and

*ξ*^{k}

*ξ*^{j}=

*ξ*^{k+j}. However, power series usually have a finite radius of convergence, which could even be zero [20,21], and therefore when the amplitude of the chaos modes exceeds a certain threshold, the series (5.7) may diverge.

^{4}

### (b) Uniformly distributed random perturbations

So far, we have been dealing with SPDEs subject to Gaussian perturbations. As we have already mentioned, similar results could also be obtained for other random perturbations, for example, uniformly distributed ones. For example, let us consider the NS equations (5.1) subject to a random initial condition ** u**(

**,**

*x**t*

_{0};

*η*), where

*η*is a uniformly distributed random variable in [−1,1]. In this setting, the solution can be expanded in a series similar to (5.4), where monic Hermite polynomials are replaced by

*rescaled Legendre polynomials*

*l*

_{k}(

*η*) (see equation (A8) in the appendix). The definition of the Wick product and the stochastic derivative can be extended to

*l*

_{k}. This yields the following analogue of formula (2.2): 5.8where are

*weighted*Wick-type products defined in (A15), whereas are Malliavin-type derivatives. Thus, the Legendre version of equations (5.2) can be written as 5.9and, apparently, it has a more complicated structure than the corresponding version for Gaussian perturbations. Nevertheless, for

*Q*=0, equations (5.9) and (5.2) are identical, and they share the same lower triangular propagator (5.6), which is therefore independent of the type of perturbation.

## 6. Numerical results

In this section, we provide numerical evidence for the theory developed above for the stochastic diffusion, Burgers and NS equations. Before doing so, however, let us briefly comment on the *computational cost* associated with the Wick–Malliavin approximation. This obviously depends on the specific problem under consideration, but there are some general facts that hold true in different situations. First of all, the Wick–Malliavin approximation relies on polynomial chaos representations. Hence, the limitations of polynomial chaos, for example, with respect to high dimensions [22–24] or discontinuities [25] are inherited in the Wick–Malliavin approach as well. In particular, the *dimensionality* of the Wick–Malliavin propagator, i.e. the total number of equations for the chaos modes, is related to the number of random input variables *M* and the Wiener chaos order *P* by the well-known formula
6.1This means that we have an exponential growth of equations with both *M* and *P*. Second, the computational complexity of the Wick–Malliavin approximation is related to the *coupling structure* of the propagator. For SPDEs with quadratic nonlinearities, this is defined by the matrices (2.9) (see also figure 1 and table 1). In particular, the structure of the Wick propagator (*Q*=0) is always *lower triangular* and *very sparse*, even for a small number of random input variables. This means that even though for a fixed Wiener chaos order, the number of equations in the propagator grows exponentially with the number of random input variables, each equation can be solved independently in sequence defined by a forward substitution algorithm. Thus, the computational cost to solve Wick-type equations is comparable with the cost of solving one equation of the propagator as many times as the total number of equations. For instance, if we consider *M*=10 random input variables and Wiener chaos order *P*=8, then we have (*K*+1)=43 758 equations. In this case, the cost of solving the Wick propagator is comparable with the cost of sampling one equation of the propagator (which is linear) 43 758 times. Note that if we use the classical probabilistic collocation method to solve the full nonlinear SPDE with the same resolution, then we would need to sample it at 9^{10} quadrature points. When the first-order Malliavin derivative is included in the stochastic approximation, the corresponding Wick–Malliavin propagator is no longer lower triangular as in the Wick setting, and we have a weak coupling between all the chaos modes. In these cases, one can consider an operator-splitting method in time, where the coupled part of the system, represented by the entries above the diagonal in the matrices of figure 1 (case *Q*=1), is treated explicitly by using the chaos modes at the previous time step. This introduces a splitting error that is however dominated by the error owing to the Wick–Malliavin truncation. Similar types of operator splitting can be considered also for higher-order Wick–Malliavin approximations.

### (a) Diffusion equation

Let us consider the diffusion problem (3.1) in a one-dimensional spatial domain [0,1] and set *f*=1. The random process *G*(*x*;*ω*) defining the random diffusivity (3.5) through its eigenvalues and eigenfunctions (*λ*_{k},*ϕ*_{k}(*x*)) is assumed to be exponentially correlated, i.e.
6.2This allows us to analytically determine the Karhunen–Loève series of *G*(*x*;*ω*) (see [26]). In particular, by setting ℓ_{c}=5 in (6.2), we have that the first three eigenfunctions include 98.4% of the total energy of the process. In the following, we consider such a case, i.e. a random diffusivity defined in terms of *M*=3 random variables. We also set the Wiener chaos expansion order to *P*=8. With these parameters, the propagator (3.8) has 165 equations and its coupling structure is shown in the last row of figure 1 for different Wick–Malliavin orders *Q*.

In figure 2, we compare the standard deviation of the solution with the full diffusion problem (3.1) with the corresponding standard deviation obtained from the Wick–Malliavin approximation (3.3). We find a very good agreement, even at very low orders of *Q*. Next, we study the convergence of the Wick–Malliavin approximation with respect to the noise amplitude *σ*. This is done in figure 2*c*, where we plot the relative difference in the *L*_{2} norm between the solution to the Wick–Malliavin approximation of order *Q* (3.3) and the full diffusion equation (3.1). Note that we obtain a fourth-order convergence rate for *Q*=1 and a sixth-order convergence rate for *Q*=2. This extends the results of Wan *et al.* [8] who demonstrated that the Wick approximation (3.4), i.e. the Wick–Malliavin approximation of order *Q*=0, has a second-order converge rate. The improvement in the convergence rate for *Q*=1 and *Q*=2 comes at the expense of a weak coupling between the equations of the propagator, which is no longer lower triangular as in the Wick setting.

### (b) Burgers equation

We first construct a benchmark stochastic solution to (4.1) by using a high-order multi-element probabilistic collocation method (ME-PCM) [27]. The ME-PCM benchmark solution allows us to study the accuracy of the Wick–Malliavin approximation (4.2) on a reliable basis, as a function of the noise amplitude *σ*. This is done in figure 3 for *M*=1 (one random variable) and different *σ*. In particular, we compare the errors in the *L*_{2} norm between the standard deviation fields obtained from the linear perturbation method (equations (4.7)), the Wick–Malliavin approximation (equation (4.2)) and the corresponding ME-PCM benchmark solution. We note that the error decreases monotonically with the order *Q* of the Wick–Malliavin approximation. In addition, the first-order perturbation method provides results that are nearly superimposed on the Wick approximation (*Q*=0). We remark that the computational cost of the Wick approximation and the linear perturbation method is comparable because both systems can be solved by forward substitution and they have exactly the same number of chaos modes. The convergence rate with respect to the noise amplitude *σ* is studied in figure 3*c*, where we plot the relative difference in the *L*_{2} norm between the solution to the Wick–Malliavin approximation of order *Q* and the full Burgers equation. Remarkably, we obtain a second-order convergence rate for *Q*=0, a fourth-order convergence rate for *Q*=1 and a sixth-order convergence rate for *Q*=2, in complete analogy with the results of the diffusion problem. In figure 4, we consider a higher-dimensional forcing term defined in terms of *M*=4 random variables and show convergence of the Wick–Malliavin approximation. We remark that for Wiener chaos order *P*=8 and *M*=4, the propagator (4.4) has 495 equations. We also note that the Wick approximation in this case is very sparse. In fact, the relative number of zero entries in the coupling matrix (2.9) is about 95% (table 1). Nevertheless, we see that, in this case, the Wick approximation can reproduce the standard deviation field within a relative error of about 10^{−3}.

### (c) Navier–Stokes equations

Next, we present fluid mechanics applications for two different prototype flow problems, i.e. two-dimensional incompressible flow past a circular cylinder at Reynolds number *Re*=100 subject to random inflow boundary condition with Gaussian distribution and a double shear layer at Reynolds number *Re*=5000 subject to a uniformly distributed perturbation in the vertical component of the initial velocity field. The computational domains and the boundary conditions of these prototype problems are sketched in figure 5.

#### (i) Flow past cylinder

Let us first consider open flow past a circular cylinder with the dimensionless inflow boundary condition assumed to be random and depending only on one normal random variable *ξ*(*ω*), i.e. we set *u*=1+*σξ*, *v*=0 where the parameter *σ*≥0 controls the amplitude of the inflow perturbation. The initial flow condition is assumed to be deterministic and coincident with a snapshot of the deterministic time-periodic wake at Reynolds number *Re*=100. The numerical discretization of the NS and the WMNS equations is performed by using the high-order spectral element method described by Karniadakis & Sherwin [28].

Our first finding is that the mean flow is not significantly affected by the truncation of the Wick–Malliavin expansion, as revealed by comparing the velocity fields obtained by the full NS equations (5.1) and the WNS equations (5.3). This is shown in figure 6, where we plot the *L*_{2} errors of the mean streamwise velocity component of the WNS equations relative to the full NS equations. This suggests that the deterministic NS equations can provide a good approximation of the mean flow in the short-term dynamics we are considering here.^{5}

On the other hand, the standard deviation is significantly different if computed from the NS or the WNS equations. This is demonstrated in figure 7, where we show one time snapshot of the standard deviation of the streamwise velocity component. We see that the standard deviation of the WNS solution grows much faster than the one of the full NS solution. In particular, there are several localized spatial regions downstream where such growth is more pronounced. This phenomenon is significantly mitigated when the first few Malliavin derivatives are included within the propagator. This can be seen in the error plots of figure 8*a*, where it is shown that Wick–Malliavin approximations of increasing order converge monotonically to the reference NS solution.

We also investigate the effect of renormalization of the WNS solution with the Catalan numbers (see §5). This is done in figures 7 and 9*a*, where we compare the error after renormalization based on an optimal (in the second-order moment sense) rescaling as well as on the Catalan numbers. The optimal rescaling is constructed by minimizing the *L*_{2} error between the variance of the solution to the full NS and the WNS equations. The effects of the renormalization are also shown in figure 9*b*, where we compare the standard deviation of the streamwise velocity component along the crossline *x*=2 (two diameters behind the cylinder) at time *t*=4.5. It is seen that the rescaling of the Wick solution by the Catalan numbers, in this case, is nearly optimal. Note that such rescaling is given *a priori*, i.e. it has not been tuned to fit the NS solution. In addition, the fact that the renormalization of the WNS solution by the Catalan numbers yields a better approximation of the NS solution is a numerical result, which we are currently investigating theoretically.

#### (ii) Double shear layer

We now consider an initial-value problem involving two fluid layers flowing in opposite directions within a periodic two-dimensional box of length two. Different from the flow past a cylinder, in this example, we set *Re*=5000 and we consider a random initial flow velocity , , *ϵ*=40, where *η* is a *uniform* random variable with a range in [−1,1] (figure 5). The perturbation in the vertical velocity component triggers the shear-layer instability, and consequently the two fluid layers roll up into two vortexes with trailing arms. The Wick–Malliavin approximation to the stochastic NS equations, in this case, is given by equations (5.9). In figure 10, we show the mean and the standard deviation of the vorticity field
6.3at time *t*=1 for *σ*=0.2. The large value of the standard deviation suggests that the small perturbations in the initial conditions have a large influence on the flow. The accuracy of the Wick–Malliavin series in representing this type of flow is studied in figure 10*c*, where we show convergence to the benchmark stochastic vorticity as we increase the order *Q*.

## 7. Summary and discussion

We presented new results regarding the accuracy and the computational efficiency of the Wick–Malliavin approximation to nonlinear SPDEs. Specifically, we studied several fundamental equations of mathematical physics, i.e. the diffusion equation with inhomogeneous random diffusivity, the stochastic Burgers driven by random noise and the NS equations subject to random boundary and random initial conditions.

We have shown that the computational complexity of the Wick–Malliavin approximation is related to the dimensionality and the coupling structure of the related propagator, i.e. the system of equations for the polynomial chaos modes of the solution. In particular, the propagator associated with the Wick approximation has a very sparse lower triangular structure that is seemingly universal, i.e. independent of the type of noise. This has significant computational advantages because each equation of the Wick propagator can then be solved independently in a sequence by a forward substitution algorithm.

On the other hand, the Wick approximation is a rather draconian stochastic truncation of nonlinearities appearing in an SPDE, and hence we addressed the question of accuracy, that is how far the Wick solution deviates from the solution to the full problem. Our numerical results for diffusion, Burgers and NS equations suggest that the Wick approximation provides a model that is relatively accurate for reasonably small noise levels. To improve the accuracy, we also introduced new higher-order stochastic approximations via Wick–Malliavin series expansions. In this generalized framework, Wick-type SPDEs can be regarded as a zero-order truncation of the Wick–Malliavin series. Higher-order approximations can be constructed in a systematic way by including Malliavin derivatives of different orders. We demonstrated through numerical simulations that the Wick–Malliavin approximation converges rapidly to the solution to the full problem as we increase the number of expansion terms. Specifically, we obtained second-, fourth- and sixth-order convergence rates in the noise amplitude as we increase the Malliavin derivative order from 0 (Wick approximation) to 2 (level-two Wick–Malliavin approximation). The improvement in accuracy and convergence rate, however, comes at the expense of a weak coupling between the equations of the propagator, which is no longer lower triangular as in the Wick setting. Both the theory and the algorithms we presented here can be readily extended to other nonlinear SPDEs and more general noises.

## Funding statement

This work was supported by OSD-MURI grant no. FA9550-09-1-0613.

## Appendix A. Wick–Malliavin rules for univariate polynomials

Here, we briefly review the Wick–Malliavin rules for univariate Hermite polynomials and also present generalizations to other families of orthogonal polynomials, in particular, Legendre polynomials.

**(a) Monic Hermite polynomials**

Let us consider the set of monic Hermite polynomials defined by the Rodrigues formula
A1By using results of linearization theory [15], we express the direct product of *h*_{i} and *h*_{j} as
A2The Wick product between two monic Hermite polynomials is defined as [7,14]
A3and, as easily seen, it represents the zero-order expansion term in (A2). The Malliavin derivatives of univariate monic Hermite polynomials can be defined as [11]
A4The Wick product of two Malliavin derivatives of the same order is
A5A comparison between (A5) and (A2) gives us the following product expansion formula in terms of Wick–Malliavin operators:
A6By using this expression, it is straightforward to show that if *X* and *Y* are two regular functions of a Gaussian variable, then we have
A7

**(b) Rescaled Legendre polynomials**

Let us define the *rescaled* Legendre polynomials
A8and [·]_{n} is the Pochhammer symbol
A9By using the results of Park & Kim [29], we obtain the direct product formula
A10The zero-order (with respect to *p*) approximation of the product *l*_{m}*l*_{n} is the polynomial *l*_{m+n}, which can be considered as the definition of a Wick product for rescaled Legendre polynomials (A8), i.e.
A11Based on this definition, we can construct the Wick calculus for rescaled Legendre polynomials. By following the same steps as in the case of monic Hermite polynomials, we arrive at
A12where and we have set
A13If *X*(*η*) and *Y* (*η*) are two functions of a uniform random variable *η*, then it follows from (A12) that
A14where the *weighted* Wick products are defined as
A15This definition gives us an interpretation of weighted Wick products for rescaled Legendre polynomials in terms of the Wick product arising from the zero-order expansion (A11). In particular, for *p*=0, we easily obtain that . Note also that for *p*>0, is *not* associative. If we apply this generalized framework to the case of Hermite polynomials of Gaussian random variables, then we obtain , for all *p*.

## Footnotes

↵1 A rigorous definition of the Wick product can be found in Janson [7, ch. 3].

↵2 The diffusion coefficient (3.5) coincides with the Wick exponential .

↵3 This happens also for the Wick–Burgers equation we discussed in §4. However, any finite truncation of the Wiener chaos expansion yields a system with finite energy.

↵4 Note that the equivalence between power series solutions and Wick propagators can be formally extended to even more general sets of basis functions in (5.7), such that

*g*_{k}()*ξ**g*_{k}()=*ξ**g*_{k+j}().*ξ*↵5 In fact, we recall that the first equation in the WNS propagator (5.6), coincides with the deterministic NS equations.

- Received January 3, 2013.
- Accepted July 15, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.