New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
A stochastic modeling methodology based on weighted Wiener chaos and Malliavin calculus

Edited by George C. Papanicolaou, Stanford University, Stanford, CA, and approved June 16, 2009 (received for review March 4, 2009)
Abstract
In many stochastic partial differential equations (SPDEs) involving random coefficients, modeling the randomness by spatial white noise may lead to illposed problems. Here we consider an elliptic problem with spatial Gaussian coefficients and present a methodology that resolves this issue. It is based on stochastic convolution implemented via generalized Malliavin operators in conjunction with weighted Wiener spaces that ensure the ellipticity condition. We present theoretical and numerical results that demonstrate the fast convergence of the method in the proper norm. Our approach is general and can be extended to other SPDEs and other types of multiplicative noise.
Elliptic partial differential equations perturbed by spatial noise provide an important stochastic model in applications, such as diffusion through heterogeneous random media (1), stationary Schrodinger equation with a stochastic potential (2), etc. A typical model of interest is the following stochastic partial differential equation (SPDE): where ω indicates randomness. Recently, this model has been actively investigated in the context of uncertainty quantification (UQ) for mathematical and computational models, see, e.g., refs. 3 and 4 and the references therein. The diffusion coefficient a(x,ω) is typically modeled by the Karhunen–Loève type expansion a(x,ω) = ā(x) + ɛ(x). Here, ā(x) is the mean and ɛ(x) = ∑_{k≥1}σ_{k}(x)ξ_{k} is the noise term, where σ_{k}(x) are deterministic matrices and ξ := {ξ_{i}}_{i≥1} is a set of uncorrelated random variables with zero mean and unit variance. For the problem of Eq. 1 to be well posed, the ellipticity condition is required. However, if the noise ɛ(x) is Gaussian, (or any other type with infinite negative part), the standard ellipticity condition would not hold, no matter how small the variance of ɛ(x,ω) is. Various modifications of the aforementioned setting have been used to mitigate this problem. For example, Gaussian models were replaced by distributions with finite range (e.g., uniform distribution). We note that from the statistical and, by implication, UQ perspectives, Gaussian perturbations are of paramount importance and should be modeled judiciously.
In this article, we propose a wellposed modification of the model of Eq. 1 with Gaussian noise ɛ = ɛ(x). More specifically, we replace Eq. 1 by the following SPDE: where “δ_{ɛ(x)}” denotes the Malliavin divergence operator (MDO) with respect to Gaussian noise ɛ(x) (see, e.g., refs. 5 and 6 and the next section). The MDO δ_{ɛ}(ζ) is a stochastic integral and specifically, in the present setting, a convolution of the integrand ζ with the driving Gaussian noise ɛ. Although replacing the ordinary product by stochastic convolution may be surprising at first, it should be appreciated that physical systems rarely produce an instant response to sharp inputs such as multiplicative forcing. Therefore, in modeling systems subject to sharp perturbations, convolutions often become a model of choice.
Historically, stochastic convolutionbased models were used to bypass the exceeding singularity of models with (literally understood) multiplicative noise. In fact, the idea of reduction to MDO could be traced back to the pioneering work of K. Itô on stochastic calculus and stochastic differential equations. Specifically, in his seminal work (7), Itô has replaced the “product” model by the “stochastic convolution” model where ∫_{0}^{t} b(u(s))dW(s) is the famous Itô's stochastic integral. In this article, we extend Itô's idea to the infinite dimensional stationary system in Eq. 1 and replace it with Eq. 2. In contrast to Itô's SDEs, elliptic equations like 1 or 2 are not causal. By this reason, we replace Itô's integral by the MDO \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $${\delta }_{\stackrel{.}{W}}$$ \end{document} . The latter, in this instance, is equivalent to anticipating (noncausal) Skorohod integral, (see ref. 6).
A general theory of bilinear elliptic SPDEs that includes, in particular, Eq. 2 was developed recently in ref. 8. In particular, it was shown that to ensure wellposedness of Eq. 2 it suffices to assume positivity of ā(x) rather than of a(x,ω) = ā(x) + ɛ(x). This approach is based on Malliavin calculus and Wiener chaos expansion (WCE) with respect to Cameron–Martin basis (see ref. 9 and next section).
The Cameron–Martin basis consists of random variables ξ_{α}, where α = (α_{1},α_{2},…) is a multiindex with nonnegative integer entries (see the next section for more detail). The WCE solution of Eq. 2 is given by the series u(x) = ∑_{α∈J}u_{α}(x)ξ_{α}, where u_{α} = 𝔼[uξ_{α}]. One can view the Cameron–Martin expansion as a Fourier expansion that separates randomness from the deterministic “backbone” of the equation. It was shown in ref. 8 that the WCE coefficients u_{α}(x) for Eq. 2 verify a lower triangular system of linear deterministic elliptic equations (see Eq. 19 below): We refer to this system as the (uncertainty) propagator. Under very general assumptions, the propagator is equivalent to SPDE 2 in the same way as the set of coefficients of a Fourier expansion is equivalent to the underlying function.
In this article, we develop a numerical method for solving Eq. 2 with highorder discretization in the physical space and weighted WCE in the probability space. We conduct a theoretical and numerical investigation of the propagation of the truncation errors and illustrate the results by numerical simulations. In particular, an a priori error estimate is presented to demonstrate the convergence of the developed numerical algorithm. We also carry out numerical and theoretical comparisons of the model of Eq. 2 with direct multiplication model of Eq. 1 for positive (lognormal) coefficient a(ω). To facilitate the theoretical comparison of these models, we have developed a generalization of the MDO δ_{ɛ} (based on Gaussian noise ɛ) to a divergence operator with respect to a class of noises including nonlinear transformations of Gaussian noise (e.g., lognormal).
In addition to its physical meaning, there are other mathematical and computational reasons for employing the “convolution” model of Eq. 2 instead of the “multiplication” model of Eq. 1, specifically: (i) The convolution model is computationally efficient due to the lower triangular (in fact, bidiagonal) structure of the propagator; its computational complexity is linear. One could also formally define an approximate WCE solution for the multiplication model of Eq. 1; however, its propagator would be a full system with quadratic computational complexity. (ii) The variance of WCE solutions for both models is infinite. However, the blowup of the convolution model is controllable, in that the WCE solution can be effectively rescaled by simple weights r_{α} such that ∑_{α}r_{α}^{2}u_{α}^{2} < ∞ (see Theorem 3). We remark that the blowup in both models is inevitable even if the perturbation ξ is 1dimensional (see example in the next section). The blowup of the multiplication model of Eq. 1 is much more severe and could hardly be controlled effectively. In fact, we are not aware of any systematic treatment of Eq. 1.
There are alternative ways to address Eq. 2 and similar bilinear equations. In particular, one could attack Eq. 2 using Hida's whitenoise infinite dimensional calculus (see refs. 10 and 11). In Hida's approach, convolution is interpreted in terms of the Wick product, an operator closely related to stochastic integrals. Convolution models for SPDEs with positive (lognormal) and normal random coefficients have been studied extensively (see, e.g., refs. 2 and 12–14) by means of whitenoise analysis. The whitenoise approach relies on builtin spaces of stochastic distributions known as Hida and Kondratiev spaces (see, e.g., refs. 15 and 16). The Malliavin calculus is more flexible, and in applications to SPDEs it allows us to build solution spaces optimal for the equation at hand (see, e.g., refs. 8 and 17). In this article, we take advantage of this important feature of Malliavin calculus to obtain more powerful numerical approximation schemes and substantially more precise estimates of the convergence rates than those reported in literature on whitenoise analysis.
Malliavin Calculus and Elliptic SPDEs
Let 𝔽 := (Ω, F, P) be a probability space, where F is the σalgebra generated by ξ := {ξ_{i}}_{i≥1} and U be a real separable Hilbert space. Given a real separable Hilbert space X, we denote by L_{2}(𝔽; X) the Hilbert space of squareintegrable Fmeasurable Xvalued random elements f. When X = ℝ, we write L_{2}(𝔽) instead of L_{2}(𝔽;ℝ). A formal series \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\stackrel{.}{W}={\sum }_{k\ge 1}{\xi }_{k}{u}_{k}$$ \end{document} , where {u_{k},k ≥ 1} is a complete orthonormal basis in U, is called Gaussian white noise on U. Specifically, for the elliptic SPDE we consider here the proper spaces are X = H^{1}(D) and U = L_{2}(D).
To explain the exact nature of solution to SPDE 2, we recall the notion of Cameron–Martin basis. Let J be the collection of multiindices α with α = (α_{1},α_{2},…) so that α_{k} ∈ {0,1,2,…} and α := ∑_{k≥1} α_{k} < ∞. For α, β ∈ J, we define α + β = (α_{1} + β_{1},α_{2} + β_{2},…), a = ∑_{k≥1}α_{k}, and α! = ∏_{k≥1} α_{k}!. By ɛ_{i} we denote the multiindex of length 1 and with the single nonzero entry at position i: i.e., the k th coordinate of ɛ_{i} is 1 if k = i and 0 if k≠i. Define the collection of random variables Ξ = {ξ_{α},α ∈ J} as follows: \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $${\xi }_{\alpha }={\prod }_{k\ge 1}({H}_{{\alpha }_{k}}({\xi }_{k})/\sqrt{{\alpha }_{k}!})$$ \end{document} , where H_{n}(x) are 1dimensional Hermite polynomials of order n.
Theorem 1. [Cameron–Martin (9)] The set Ξ is an orthonormal basis in L_{2}(𝔽): if η ∈ L_{2}(𝔽) and η_{α} = 𝔼[ηξ_{α}], then η = ∑_{α ∈ J} η_{α}ξ_{α} = ∑_{α ∈ J}η_{α}H_{α}(ξ)/ \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\sqrt{\alpha !}$$ \end{document} and 𝔼[η^{2}] = ∑_{α ∈ J}η_{α}^{2}.
Next, we introduce the Malliavin derivative (MD) and the MDO. To make these notions more transparent, we begin with the simplest setting: we will differentiate and integrate ξ_{α} with respect to a single normal random variable, e.g. ξ_{k}, the kth coordinate of ξ. The MD D_{ξk} and divergence operator δ_{ξk} are defined, respectively, by
In the literature on quantum physics (see, e.g., ref. 18) the operators δ_{ξk} and D_{ξk} are often called creation and annihilation operators, respectively. As intuition suggests, the MDO of D_{ξk}(ξ_{α}) recovers the integrand ξ_{α}. Indeed, for α_{k} > 0, α_{k}^{−1}δ_{ξk}(D_{ξk}(ξ_{α})) = ξ_{α}. Also, MD and MDO can be extended to L_{2}(𝔽;X). For example, for f ∈ L_{2}(𝔽;X⊗U), \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $${\delta }_{\stackrel{.}{W}}(f)$$ \end{document} is defined as the unique element of L_{2}(𝔽;X) with the property 𝔼[φ \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $${\delta }_{\stackrel{.}{W}}(f)$$ \end{document} ] = 𝔼[(f,D_{ \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\stackrel{.}{W}$$ \end{document} }φ)_{U}] for every φ such that φ ∈ L_{2}(𝔽) and D_{ \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\stackrel{.}{W}$$ \end{document} }φ ∈ L_{2}(𝔽;U) (see, e.g., ref. 17).
Before proceeding with the SPDE 2, let us consider a simple example to shed some light on the structure of the spaces within which we could expect existence of solutions of this equation.
Example 1. Let u be a solution of equation where ξ ∼ N(0,1). Simple calculations show that ∥u∥_{L2(𝔽)}^{2} = ∑_{k = 1}^{∞} u_{n}^{1} = ∞, where u_{n} = ∑_{n≥0}𝔼[uH_{n}(ξ)]/ \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\sqrt{n!}$$ \end{document} . On the other hand, taking a weighted norm with weights r_{n} = (n!)^{−1/2}2^{−qn/2},q > 0, we get
The above example demonstrates that even simple stationary equations do not have solutions with finite second moments. To overcome this obstacle one should introduce an appropriate weighted version of the solution space. Clearly, introduction of weights amounts to rescaling of the stochastic Fourier (Wiener chaos) representation of the solution. Below, we discuss briefly the construction of weighted spaces.
Let R be a bounded linear operator on L_{2}(𝔽) defined by Rξ_{α} = r_{α}ξ_{α} for every α ∈ J, where the weights {r_{α}, α ∈ J} are positive numbers. In what follows, we will identify the operator R with the corresponding collection {r_{α}, α ∈ J}. The inverse operator R^{−1} is defined as R^{−1}ξ_{α} = r_{α}^{−1}ξ_{α}. The elements of RL_{2}(𝔽; X) can be identified with a formal series ∑_{α ∈J}f_{α}ξ_{α}, where ∑_{α∈J}∥f_{α}∥_{X}^{2}r_{α}^{2} < ∞. Clearly, RL_{2}(𝔽; X) is a Hilbert space with respect to the norm f_{RL2(𝔽; X)}^{2} := Rf_{L2(𝔽; X)}^{2}. We define the space R^{−1}L_{2}(𝔽; X) as the dual of RL_{2}(𝔽; X) relative to the inner product in the space L_{2}(ℝ; X). For f ∈ RL_{2}(𝔽; X) and g ∈ R^{−1}L_{2}(𝔽) we define the scalar product where 〈〈f,g〉〉 = ∑_{α ∈ J}f_{α}g_{α}, with g = ∑_{α ∈ J}g_{α}ξ_{α}.
Remark 1. One could readily check that u = 1 + ξ · u, i.e., the multiplication version of Eq. 6, is much more complicated and blows up much faster than Eq. 6.
To address SPDEs driven by lognormal and other important types of random perturbations, we need to develop a bilinear symmetric version of MDO with white noise \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\stackrel{.}{W}$$ \end{document} replaced by an arbitrary nonlinear transformation of Gaussian random variables. To begin with, assume that vector ξ consists of a single Gaussian random variable ξ ∼ N(0, 1). Let ξ_{m} = H_{m}(ξ)/ \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\sqrt{m!}$$ \end{document} . Define a ntuple MDO by induction: δ_{ξ}^{⊗1}(ξ_{m}) := δ_{ξ}(ξ_{m}) and δ_{ξ}^{⊗n}(ξ_{m}) := δ_{ξ}δ_{ξ}^{⊗(n−1)}(ξ_{m}) for n > 1. Let δ_{ξn}(ξ_{m}) := δ_{ξ}^{⊗n}(ξ_{m})/ \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\sqrt{n!}$$ \end{document} . It is readily checked that Now, for ξ := (ξ_{1}, ξ_{2}, …) and multiindices α and β, define \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $${\delta }_{{\xi }_{\alpha }}({\xi }_{\beta }):={\prod }_{k=1}^{\infty }{\delta }_{{\xi }_{k}}^{\otimes {\alpha }_{k}}({\xi }_{{\beta }_{k}}({\xi }_{k}))$$ \end{document} . The formula of Eq. 9 translates into the multidimensional case as follows:
Remark 2. (Wick product) The Wick product (◇), can be defined as follows: for α, β ∈ 𝕁, H_{α}(ξ) ◇ H_{β}(ξ) := H_{α+β}(ξ), where H_{γ}(ξ) := ∏_{k≥1} H_{γk}(ξ_{k}). It follows from Eq. 10 that H_{α}(ξ) ◇ H_{β}(ξ) = δ_{Hα(ξ)}(H_{β}(ξ)).
Theorem 2. If θ, η ∈ RL_{2}(𝔽; X), set δ_{ξα}(θ) := ∑_{β∈J} θ_{β}δ_{ξα} (ξ_{β}) and δ_{η}(θ) := ∑_{α∈J} η_{α}δ_{ξα} (θ). Then δ_{η} is a bounded linear operator on R′L_{2}(𝔽; X) and
The proof follows from Eq. 10.
Example 2. Let ξ ∼ N(0, 1) and c is a constant. Let \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\theta =exp(c\xi \frac{{c}^{2}}{2})$$ \end{document} and \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\eta =exp(c\xi \frac{{c}^{2}}{2})$$ \end{document} . Then, δ_{η}(θ) = 1. This relation is important for the comparison of the product of Eq. 1 and convolution of Eq. 2 with lognormal coefficient a(ω) (see Numerical Results).
Elliptic SPDE Model.
Let us rewrite Eq. 2 as an SPDE in the sense of Malliavin calculus: where D denotes the physical domain, \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\stackrel{.}{W}$$ \end{document} is white noise on U = L_{2}(D), the Hilbert space of square summable sequences of real numbers, and Mu = ∑_{k≥1}M_{k}u ⊗ u_{k}. We note that u_{k} in L_{2}(D) can be identified with ɛ_{k}, see Eq. 5. For simplicity, we assume that g = 0 and that f is deterministic. We can now rewrite the term \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $${\delta }_{\stackrel{.}{W}}(Mu)$$ \end{document} as Everywhere below we assume that: (i) the functions a_{ij}(x) and σ_{ij}^{k}(x) are measurable and bounded in the closure D̄ of D; and (ii) there exist positive numbers A_{1}, A_{2} such that A_{1}y^{2} ≤ a_{ij}(x)y_{i}y_{j} ≤ A_{2}y^{2} for all \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$x\in \overline{D}$$ \end{document} and y ∈ ℝ^{d}.
Definition 1. A solution to Eq. 12 is a random element u ∈ H_{0}^{1}(D) such that the equality where f ∈ H^{−1}(D) and
Analytical issues related to Eq. 15 have been investigated in ref. 8. In particular, the following result holds:
Theorem 3. There exists an operator R and a unique solution u ∈ RL_{2}(𝔽; H_{0}^{1}(D)) of Eq. 12 such that:

(i) The operator R is defined by the weights r_{α} given by where the numbers q_{k} are chosen so that the “renormalization condition” holds, and C_{k} are defined by

(ii) With the definition of r_{α} given in (i), u_{α} satisfy where the constant C_{A} is defined by
Because the expectation of the MDO is zero, we have which is the unperturbed (deterministic) version of elliptic Eq. 1. By substituting the WCE u = ∑_{α∈J}u_{α}ξ_{α} into Eq. 15 and performing a Galerkin projection in the probability space, we can establish the equivalence between Eqs. 2 and 15 and derive the following uncertainty propagator: which is a system of deterministic partial differential equations (PDEs). In the next section, we will develop an efficient numerical algorithm to solve Eq. 19.
Numerical Method
We will employ WCE in the probability space and a highorder finite element discretization in the physical space. Our method exploits the recursive structure of the propagator of Eq. 19. This remarkable property of the propagator significantly reduces the complexity of the algorithm and provides opportunities for improving numerical efficiency. For example, the coefficients u_{α}, satisfying α = p, can be solved in parallel because they all rely only on already computed coefficients u_{β} with β < α and α − β = 1.
For simplicity and without loss of generality, we consider a 2dimensional physical domain D. Let T_{h} be a family of triangulations of D with straight edges and h the maximum size of elements in T_{h}. We assume that the family is regular, in other words, the minimal angle of all the triangles is bounded from below by a positive constant. We define the finite element space as and where F_{K} is the mapping function for the element K which maps the reference element R to element K and P^{ \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\widehat{p}$$ \end{document} }(R) denotes the set of polynomials of degree up to \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\widehat{p}$$ \end{document} over R. We assume that v_{h}_{ΓD} = 0, ∀v_{h} ∈ V_{h}. Thus, V_{h} is an approximation of H_{0}^{1}(D) based on piecewise polynomials. There exist many choices of basis functions on the reference elements, such as htype finite elements (19), spectral/hp elements (20, 21), etc.
Let J_{M, p} be a finite dimensional subset of J given by J_{M, p} := {αα ∈ ℕ_{0}^{M}, α ≤ p, p ∈ ℕ} and R_{M ,p} be a given set of weights r_{α}, α ∈ J_{M, p}. We define
Then, the stochastic finite element method (sFEM) can be defined as: Find u_{h}^{M, p} ∈ V_{h} ⊗ V_{c} such that for any v ∈ V_{h} ⊗ V_{c}^{−1} where
Note that we truncate the expansion of the white noise up to M terms and the WCE up to polynomial order p.
Let H = H_{0}^{1}(D), and \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $${r}_{\alpha }={q}^{\alpha }/\sqrt{\left\alpha \right!}$$ \end{document} as in Theorem 3. The main result regarding convergence of the stochastic finite element method (sFEM) is given by the following theorem:
Theorem 4. For u ∈ RL_{2}(𝔽; H) ∩ RL_{2}(𝔽; H^{m+1}(D)), the error of approximation of the sFEM is given by where the constant C is independent of h, \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\widehat{q}$$ \end{document} = ∑_{k≥1}C_{k}^{2}q_{k}^{2} < 1, \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\widehat{q}$$ \end{document} _{W} ∑_{k>M}C_{k}^{2}q_{k}^{2}, and C_{k} are constants defined in Theorem 3.
We remark that the 3 main components of the error O(h^{m}), O( \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\widehat{q}$$ \end{document} _{W}), and O( \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\widehat{q}$$ \end{document} ^{p+1}) are due to the finite element discretization, the approximation of white noise, and truncation of the WCE, respectively. We also note that spectral convergence is obtained in the weighted norm  · _{RL2(𝔽; H)}. Next, we present a sketch of the proof whereas all technical details can be found in supporting information (SI) Appendix.
The approximation error can be decomposed as where u_{α} and \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\widehat{u}$$ \end{document} _{α} are the coefficients of chaos expansions of u and u_{h}^{M, p}, respectively. Correspondingly, we obtain To obtain the error contribution I_{1}, we need to estimate the finite element approximation error u_{α} − \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\widehat{u}$$ \end{document} _{α}. Since u_{α} depends on u_{β} with α − β = 1, we should consider the error propagation in the uncertainty propagator. The key observation is that for α > 0, the following equations are satisfied in the weak sense from which we can obtain a recursive inequality for the error propagation as We then use results from the approximation theory to bound the finite element approximation error. The error contribution from I_{2} is from the truncation of WCE and the approximation of white noise, which can be estimated using Theorem 3. More details are given in SI Appendix.
Numerical Results
We consider the 2dimensional stochastic elliptic problem where 𝔼[a](x) denotes the mean field of coefficient and \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\stackrel{.}{W}(x)$$ \end{document} the random perturbation. For simplicity, we choose the physical domain D = (0,1)^{2}, 𝔼[a](x) = 1, and f(x) = 1.
To represent the white noise on L_{2}(D), we select the following orthonormal basis Hence, we approximate the white noise as For convenience, we use an 1dimensional index and rewrite the above equation as Let \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$u(x)={\sum }_{\left\alpha \right=0}^{p}{u}_{\alpha }{\xi }_{\alpha }$$ \end{document} be the truncated WCE of the solution up to polynomial order p. The uncertainty propagator takes the form where f_{α} = 0 if α > 0, because f(x) is assumed to be deterministic, and the operator M_{k} is defined as M_{k}u =−∇ · (w_{k}(x)∇u). To choose proper r_{α} for the weighted Wiener chaos space, we need to specify the constant C_{k} defined in Theorem 3. For our case, we have Here, w_{k}(x)_{L∝(D)} is uniformly bounded and A_{1} = 1 because 𝔼[a](x) = 1. Hence, the weights can be defined as \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $${r}_{\alpha }={q}^{\alpha }/\sqrt{\left\alpha \right!}$$ \end{document} with \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $${q}^{\alpha }={\prod }_{k=1}^{\infty }{q}_{k}^{{\alpha }_{k}}$$ \end{document} , if q_{k} is chosen as q_{k} = 1/(k + 1)C_{k}.
We now examine the convergence of the sFEM. We employ the spectral/hp element method to solve the uncertainty propagator. The physical domain (0,1)^{2} is discretized in 64 quadrilateral elements, where 12thorder piecewise polynomials are used in each element. Numerical tests show that errors associated with the physical discretization are close to machine accuracy and hence negligible. If we fix the number of random variables in the approximation of white noise and increase the polynomial order of WCE, the dominant error in the error estimate of Eq. 21 should be the truncation error of the WCE. In Fig. 1, we plot the convergence behavior of the approximate white noise in the weighted norm RL_{2}(𝔽,L_{2}(D)) with respect to M. We also plot the errors of the approximate solution in the weighted norm RL_{2}(𝔽,H_{0}^{1}(D)) with respect to the polynomial order of WCE. For the numerical simulation, we choose M = 21, corresponding to m + n ≤ 5. The errors are approximated as u_{h}^{M, p+1} − u_{h}^{M, p}_{RL2(𝔽; H01(D))}. We see that spectral convergence is obtained, which is consistent with the error estimate of Eq. 21.
Model Comparison.
Next, we present a comparison of the convolution and multiplication models by considering the following 1dimensional problem where “*” indicates stochastic convolution or ordinary multiplication. We choose f(x) = 1, and we consider first the simple model \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $${K}_{I}(x,\xi )=exp(c\xi \frac{1}{2}{c}^{2})$$ \end{document} , where ξ ∼ N(0,1) is a normal random variable and c is a constant indicating the degree of perturbation. In particular, K_{I}(x) can be regarded as a simplified version of the lognormal noise \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $${\hbox{ e }}^{\stackrel{.}{W}}$$ \end{document} . In addition, we consider a second model K_{II} with spacedependent lognormal noise of the type where ξ_{i}, i = 1,2,3, are normal random variables. In other words, we take the first 3 modes of white noise \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\stackrel{.}{W}(x)={\xi }_{1}+{\sum }_{i=2}^{\infty }\sqrt{2}cos((i1)\pi x){\xi }_{i}$$ \end{document} on L_{2}(0,1). It is easy to show that 𝔼[K_{II}(x, ξ)] = 1, and that Var(K_{II}(x,ξ)) = exp(c^{2} + 2c^{2}cos^{2}(πx) + 2c^{2}cos^{2}(2πx)) − 1. In Fig. 2, we plot the variances of solutions corresponding to 2 different models obtained from solving both the convolution and the multiplication cases; results for both K_{I} and K_{II} are shown. We observe that when the noise is small, the variances given by the 2 models are about the same; however, when the noise is large a large difference exists, with the variance given by the multiplication model increasing much faster than that given by the convolution model. For K_{I}, when its standard deviation increases from 0.1003 to 22.7379, the solution variance increases from O(10^{−5}) to O(10^{6}) for the multiplication model, and from O(10^{−5}) to O(10^{1}) for the convolution model—a 5order difference in magnitude! For K_{II} the variace achieves larger values as we include more terms in the expansion for the white noise. More results and analysis for both models are included in SI Appendix.
This exponential increase of the variance with respect to perturbation can be explained by referring to the aforementioned Example 2. Specifically, we consider the simple problem (δ_{a}(u_{x}(x,ξ)))_{x} = f(x),x ∈ (x_{0},x_{1}),u(x_{0}) = u(x_{1}) = 0 with lognormal but spaceindependent coefficient \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$a(\omega )=exp(c\xi \frac{{c}^{2}}{2})$$ \end{document} ; the corresponding solution is u(x) = θδ^{−1}f(x) with δ being the Dirichlet Laplacian on (x_{0},x_{1}) and δ_{a}(θ) = 1. Example 2 shows that \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\theta =exp(c\xi \frac{{c}^{2}}{2})$$ \end{document} . On the other hand, the solution to the corresponding multiplication model (a · v_{x}(x,ξ))_{x} = f(x) is v(x,ξ) = a^{−1}δ^{−1}f(x). Hence, the rapid increase in the variance observed in the computations is related to the ratio a^{−1}/θ, which is equal to exp(c^{2}). Clearly, v(x) = exp(c^{2})u(x) and 𝔼[v(x)^{k}] = exp(kc^{2})𝔼[u(x)^{k}]. This establishes that solutions of the above equation with lognormal permeability corresponding to multiplication model are exceedingly singular as compared to a similar equation with stochastic convolution based on MDO.
Summary and Discussion
In this article, we have developed an efficient numerical method for solving secondorder elliptic SPDEs with Gaussian coefficients. We introduced the concept of “stochastic convolution model” and employed Malliavin calculus to implement it. The stochastic solution was represented by a weighted WCE in the appropriate norm. It was shown that the coefficients of the WCE can be obtained by solving a lower triangular system of deterministic elliptic PDEs, i.e., the uncertainty propagator. We derived a a priori error estimate to measure the rate of convergence with respect to appropriately weighted L_{2} norm. We have also carried out numerical and theoretical comparisons of the model of Eq. 2 with the direct multiplication model of Eq. 1 for positive (lognormal) coefficient a(ω). To facilitate the theoretical comparison of these models, we have developed a generalization of the MDO δ_{ɛ} (based on Gaussian noise ɛ) to a divergence operator with respect to a class of noises including nonlinear transformations of Gaussian noise (e.g., lognormal).
It might be instructive to reexamine the differences and similarities between the convolution model \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$Au(x)={\delta }_{\stackrel{.}{W}}(Mu(x))$$ \end{document} and its multiplication counterpart Au(x) = Mu(x) · \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\stackrel{.}{W}$$ \end{document} (x). It is well understood that the multiplication models are exceedingly singular. We are not aware of any systematic theoretical or numerical efforts to investigate such SPDEs. It appears though that the general methodology developed in this article could be extended to address elliptic equations of this type as well. However, in the multiplication setting, the propagator is not lower triangular, and the solution spaces are expected to be much larger than in the setting of this article. In contrast to multiplication models the convolution models are much more manageable analytically and numerically and have solid physical meaning. Both models become much easier if one replaces the white noise forcing by lognormal (which is a positive noise). We have shown that in the lognormal setting the SPDE has reasonable solutions for both models. Still, as the variance of the noise grows, the variance of the solution of the multiplication equation scales up exponentially as compared with the convolution model.
In spite of the aforementioned differences, the 2 models are closely related. In fact, the convolution models could be viewed as the highest stochastic order approximations to multiplication models. Indeed, by Remark 2, δ_{Hα(ξ)}(H_{β}(ξ)) = H_{α+β}(ξ) while H_{α}(ξ) · H_{β}(ξ) = H_{α+β}(ξ) + R_{α,β}, where R_{α,β} is a linear combination of Hermite polynomials of orders lesser than α+β.
Finally, we note that the mean in the multiplication model for the problems considered here deviate greatly from the corresponding deterministic solution (see SI Appendix). Clearly, linear systems with additive noise are unbiased perturbations of their deterministic counterparts with the statistical average of a solution to randomized system coinciding with the solution of the original unperturbed system. The convolution bilinear equations considered in this article also enjoy this important property of linear systems. However, we stress that this property does not hold and should not be expected from SPDEs with nonlinear operator A or/and multiplication in the stochastic terms.
Acknowledgments
We acknowledge the helpful suggestions by the referees and the useful discussions with our students Z. Q. Zhang and C.Y. Lee. This work was supported by the National Science Foundation (NSF)/Division of Mathematical Sciences and Army Research Office grants (to B.R.) and NSF/Applied Mathematics Crosscutting–Stochastic Systems and Office of Naval Research (G.E.K.).
Footnotes
 ^{1}To whom correspondence should be addressed. Email: gk{at}dam.brown.edu

Author contributions: X.W., B.R., and G.E.K. performed research and wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.
References
 ↵
 Keller JB,
 McLaughlin D,
 Papanicolaou G
 Papanicolaou G
 ↵
 Holden H,
 Oksendal B,
 Uboe Zhang T
 ↵
 ↵
 ↵
 Malliavin P
 ↵
 Nualart D
 ↵
 Itô K
 ↵
 Lototsky S,
 Rozovskii B
 ↵
 ↵
 Hida T,
 Juo H,
 Potthoff J,
 Streit L
 ↵
 Kuo HH
 ↵
 Manouzi H,
 Theting T
 ↵
 Theting T
 ↵
 Våge G
 ↵
 ↵
 ↵
 ↵
 Glimm J,
 Jaffe A
 ↵
 Ciarlet P
 ↵
 Karniadakis GE,
 Sherwin SJ
 ↵
 Schwab C
Citation Manager Formats
Sign up for Article Alerts
Jump to section
You May Also be Interested in
More Articles of This Classification
Physical Sciences
Related Content
 No related articles found.
Cited by...
 No citing articles found.