University of Maryland Department of Computer Science TR-4786
University of Maryland Institute for Advanced Computer Studies TR-2006-10
Solving the Stochastic Steady-State Di usion Problem using
Multigrid1
Howard Elman2 and Darran Furnival3
February 20, 2006
ABSTRACT. We study multigrid for solving the stochastic steady-state di usion problem.
We operate under the mild assumption that the di usion coe cient takes the form of
a  nite Karhunen-Lo eve expansion. The problem is discretized using a  nite element
methodology using the polynomial chaos method to discretize the stochastic part of the
problem. We apply a multigrid algorithm to the stochastic problem in which the spatial
discretization is varied from grid to grid while the stochastic discretization is held constant.
We then show, theoretically and experimentally, that the convergence rate is independent
of the spatial discretization, as in the deterministic case.
1. Introduction
Mathematical models often contain partial di erential equations (PDEs). The constituent parts of
the PDE, i.e. the di erential coe cients and the source function, are most easily modeled as functions
of the spatial domain. However, uncertainty might exist as to the most appropriate functions to use
in the model. A more sophisticated model might therefore represent the di erential coe cients and
source function not only as functions on the spatial domain but also as functions on some sample
space, i.e. as random  elds. This gives rise to stochastic partial di erential equations (SPDEs).
In this paper we consider the stochastic steady-state di usion equation along with homogeneous
Dirichlet boundary value conditions. We are interested in the case when the di usion coe cient is
stochastic and the source function is deterministic, i.e. the di usion coe cient is a random  eld and
the source function is de ned on the spatial domain only. However, we also treat the source function
as a random  eld as this is required for purposes of analysis and incorporates, as a special case, the
fact that the source function may be deterministic.
We will assume the di usion coe cient to be of the form of a  nite Karhunen-Lo eve expansion.
This is common in the literature, e.g. see Babu ska, Tempone & Zouraris (2004), Ghanem & Spanos
(1991), and Xiu & Karniadakis (2002). We require the di usion coe cient to be of this form for
analytic as well as computational purposes.
We are interested in using a  nite element methodology to  nd an approximate solution to the
problem. We therefore obtain a weak formulation to the boundary value problem and proceed to look
in a  nite-dimensional subspace of the in nite-dimensional space that contains the weak solution in
order to obtain a matrix problem.
1This work was supported in part by the U.S. National Science Foundation under grant DMS0208015 and by the
the U.S. Department of Energy under grant DOEG0204ER25619.
2Department of Computer Science and Institute for Advanced Computer Studies, University of Maryland, College
Park, MD 20742. elman@cs.umd.edu.
3Applied Mathematics and Scienti c Computing Program, University of Maryland, College Park, MD 20742.
furnival@math.umd.edu
1
2
The  nite-dimensional subspace in which we look for the approximate solution will be a tensor
product of a space of functions de ned on the spatial domain and a space of functions de ned on
the sample space. For the  nite-dimensional space of functions on the spatial domain we will choose
the set of piecewise linear polynomials de ned on a triangulation. For the  nite-dimensional space
of functions de ned on the sample space we will choose the set of polynomials that are of degree no
greater than a chosen parameter.
For computational purposes the basis of the space of functions de ned on the sample space will
correspond to that obtained by employing the polynomial chaos method as developed by Ghanem
& Spanos (1991) and generalized by Xiu & Karniadakis (2002). Indeed, this method provided the
foundation on which this work was developed. However, this is not the only basis which is possible.
Another choice would be the doubly orthogonal polynomials as discussed in Babu ska et al. (2004).
Although we will not compute with the doubly orthogonal polynomials they will be vital to our
analytic argument.
Theoretically, we apply a two-grid correction scheme to solve the matrix problem. In this scheme
the spatial discretization is varied from grid to grid, giving a  ne grid on the upper level and a
coarse grid on the lower level, while the stochastic discetization is kept constant. We show that
the convergence factor of this method is independent of the spatial mesh parameter. An induction
argument can then be used to show that the multigrid algorithm obtained by applying the two-grid
correction scheme recursively has a convergence factor that is also independent of the spatial mesh
parameter. We do not give this inductive argument but note that the reasoning would be the same
as that for the analogous deterministic problem which is discussed e.g. in Braess (2001) and Elman,
Silvester & Wathen (2005).
Experimentally, we consider two problems. These are obtained by de ning the di usion coe cient
to be a Karhunen-Lo eve expansion consisting of random variables that are for the  rst problem
uniformly distributed and for the second problem normally distributed. In both these problems the
source function is set to unity. We then apply multigrid to a selection of matrix problems associated
with di erent discretization parameters, both spatial and stochastic, and tabulate the number of
iterates required for convergence. It will be seen from these tables that the number of iterates it
requires multigrid to converge remains, to within a close approximation, constant as the spatial
parameter is varied. This provides experimental evidence that the convergence factor associated
with the applied multigrid algorithm is indeed independent of the spatial mesh parameter.
2. The Stochastic Steady-State Diffusion Problem
In this section we introduce the stochastic steady-state di usion problem along with its weak
formulation and  nite element discretization. We then obtain a number of properties associated with
the discretized system. These properties are analogous to the properties of the system of equations
resulting from the deterministic steady-state di usion problem and are proven similarly. Finally,
we introduce the doubly orthogonal polynomials which can be used as a basis for the stochastic
part of the problem. The doubly orthogonal polynomials will be vital to our theoretical arguments
concerning the convergence of multigrid in x3.
2.1. Boundary Value Problem. The stochastic steady-state di usion equation with homogeneous
Dirichlet boundary value conditions is given by
  r (cru) = f in D  ;
u = 0 on @D  ;(1)
3
where D is the spatial domain,  is a sample space, c: D  !R is the di usion coe cient, and
f: D  !R is the source function. The sample space in turn belongs to a probability space
( ;F;P) where F is a  -algebra and P is a probability measure.4
The spatial domain, D, is assumed to be a two-dimensional simply connected bounded open set
with piecewise smooth boundary. In particular we take D to be the interior of a polygon.
We will let  =  m where  = (a;b). We will assume the di usion coe cient to be of the form
c(x;!) = c0(x) +
mX
r=1
p
 rcr(x) r(!)(2)
where  r:  !R are identically distributed independent random variables with zero mean and for
! = (!1;:::;!m) 2  ,  r(!) = !r. Note that the distribution of  = ( 1;:::; m) will dictate the
probability measure to be used. For example, if  r is uniformly distributed on ( 1;1) then P will
be the probability measure associated with an m-dimensional uniform distribution.
We necessarily expect the solution to be a random  eld, u: D  !R, such that for each value
of !2 the resulting PDE is satis ed in the classical sense.
We note that this problem is extensively discussed from a modeling perspective in Ghanem &
Spanos (1991) and from an analytic perspective in Babu ska et al. (2004).
2.2. Weak Formulation. In stating the weak formulation of (1) we will use tensor products of
Hilbert spaces which are de ned and discussed in Babu ska et al. (2004) and Treves (1967). Let
c 2 L1(D) L1( ) and f 2 L2(D) L2( ). The weak formulation of (1) is given by:  nd
u2H10 (D) L2( ) such that
a(u;v) = l(v) 8v2H10 (D) L2( )(3)
where5
a(u;v) =
Z
 
Z
D
cru rv;(4)
l(v) =
Z
 
Z
D
fv:(5)
The Lax-Milgram lemma can be used to show that there exists a unique solution to this problem
providing that there exist positive constants  and  such that
  c(x;!)  P-a.e.8x2D:(6)
2.3. Finite Element Formulation. We are interested in applying a  nite element methodology
to  nd an approximation to the solution of the variational problem given in x2.2. This entails a
discretization of both the spatial and stochastic parts of the problem.
The spatial domain is discretized using a triangulationT =f41;:::;4Kg. Denoting the longest
side of the t-th triangle in the triangulation as ht, we de ne the mesh parameters h = maxht and
h = minht. We also denote the smallest angle in the t-th triangle as  t. We assume that any
triangulation used belongs to a family of triangulations that is quasi-uniform, i.e. there exists  > 0
4Note that the nabla operator only operates on the spatial components of the function, that is to say, that if
f: D  !R and x = (x1;x2)2D, !2 , then
rf(x;!) =
 @f(x;!)
@x1 ;
@f(x;!)
@x2
 
:
5Note that the integral over  is with respect to the probability measure, i.e.
Z
 
Z
D
=
Z
 
Z
D
dxdP:
4
such that h> h for every triangulation in the family, and shape regular, i.e. there exists   > 0 such
that    min t for every triangulation in the family. The  nite-dimensional subspace of H10 (D) is
then taken to be S = spanf 1;:::; Ng where  k: D!R, k = 1;:::;N, are the usual piecewise
linear basis functions de ned at the nodes of T. (Here N is the number of internal nodes in the
triangulation.)
For the subspace of L2( ) we take the space of all v2L2( ) such that v is a polynomial of degree
at most p. We denote this space by T. For the purposes of analysis any basis of T is su cient. We
will use the notation T = spanf 1;:::; Mg where M = M(m;p) and is in fact given by
M = (m+p)!m!p! :(7)
The basis with which we choose to compute in x4 is that derived from the generalized polynomial
chaos method as discussed in Xiu & Karniadakis (2002). In this method the functions chosen to
be a basis for T are those polynomials from the Askey scheme of hypergeometric polynomials that
satisfy Z
 
 k l = dk kl:(8)
For example, if  r, r = 1;:::;m, are uniformly distributed on ( 1;1) then the basis of T would be
the set of m-dimensional Legendre polynomials of degree at most p.
We thus have S T  H10 (D)  L2( ) which leads to the  nite element formulation:  nd
uhp2S T such that
a(uhp;v) = l(v) 8v2S T(9)
where a( ; ) and l( ) are as in (4) and (5). This will possess a unique solution under the same
conditions that apply to the weak formulation.
2.4. Matrix Formulation. By substituting the expansion
uhp =
NX
j=1
MX
l=1
ujl j l(10)
into (9) and varying v over the basis functions of S T we  nd that we can obtain the  nite element
approximation by solving the matrix problem:  nd u2RMN such that
Au = f(11)
where
A =
2
64
A11    A1M
... ...
AM1    AMM
3
75; [A
kl]ij =
Z
 
Z
D
cr i r j k l;(12)
and
f =
2
64
f1
...
fM
3
75; [f
k]i =
Z
 
Z
D
f i k;(13)
with u = [u11;:::;uN1;:::;u1M;:::;uNM]T.
Once we have computed u then we have uhp. From this we can calculate such things as the
mean, variance, and covariance of the approximation. Formulas for these quantities in terms of the
coe cients of (10) are readily obtainable once a basis for T has been chosen.
5
2.5. Matrix and Right Hand Side Properties. We here establish some results concerning the
system matrix, A, and the right hand side vector, f, that will be required for the analysis of multigrid
in x3.
In the following, E refers to the stochastic mass matrix and B refers to the (deterministic) mass
matrix, which are de ned by
[E]kl =
Z
 
 k l; [B]ij =
Z
D
 i j;(14)
respectively. By E B we will mean the matrix Kronecker product as given by
E B =
2
64
e11B    e1MB
... ...
eM1B    eMMB
3
75:(15)
We now introduce some notation for the coe cient vector of a function in S T. Let v2S T,
then we have the expansion
v =
NX
i=1
MX
k=1
vik i k:(16)
We de ne the coe cient vector v2RMN of v by
v =
2
64
v1
...
vM
3
75; [v
k]i = vik:(17)
Theorem 1. Let f2S T with coe cient vector ^f, then f = (E B)^f, where f is as in (13).
Proof. This follows upon substituting the expansion of f into (13).  
Theorem 2. Let f2S T with coe cient vector ^f, then jjfjj2L2(D) L2( ) = ((E B)^f;^f).
Proof. This follows upon substituting the expansion of f into jj jjL2(D) L2( ).  
Theorem 3. The inequality
C1h2  ((E B)v;v)(v;v)  C2h2(18)
holds for all v2RMN, where C1 = C1(p) and C2 = C2(p).
Proof. Let v2RMN be the coe cient vector of some v2S T. Then,
((E B)v;v) =
NX
i=1
MX
k=1
[vk]i
NX
j=1
MX
l=1
[E]kl[B]ij[vl]j(19)
=
NX
i=1
MX
k=1
vik
NX
j=1
MX
l=1
vjl
Z
 
 k l
Z
D
 i j
=
KX
t=1
NX
i=1
MX
k=1
NX
j=1
MX
l=1
vikvjl
Z
 
 k l
Z
4t
 i j:
6
Let Bt be the local mass matrix and vt2R3M be a vector that contains only those values of v that
are associated with the t-th triangle in the mesh, viz.,
vt =
2
66
4
v(t)1
...
v(t)M
3
77
5; [v
(t)
k ]i = v
(t)
ik(20)
where v(t)k 2R3 contains the values of v at those vertices of the t-th triangle associated with the
k-th basis function of T. Let B be the mass matrix on the canonical triangle (with vertices (0;0),
(1;0), (0;1)). Then
((E B)v;v) =
KX
t=1
3X
i=1
MX
k=1
3X
j=1
MX
l=1
v(t)ik v(t)jl
Z
 
 k l
Z
4t
 (t)i  (t)j(21)
=
KX
t=1
2j4tj
3X
i=1
MX
k=1
[v(t)k ]i
3X
j=1
MX
l=1
[E]kl[B ]ij[v(t)l ]j
=
KX
t=1
2j4tj((E B )vt;vt)
where j4tj is the area of the t-th triangle.
For any triangle with longest side ht and smallest angle  t we have the inequality
1
4h
2
t sin t j4tj 
1
2h
2
t sin t;(22)
and given that E B is a symmetric positive de nite matrix there exists  =  (p) and  =  (p)
such that
 (vt;vt) ((E B )vt;vt)  (vt;vt)(23)
for all vt.
Using (22) and (23) along with the assumptions of quasi-uniformity and shape regularity, and
continuing with (21) we have
(E Bv;v) 
KX
t=1
1
2 h
2
t sin t(vt;vt)(24)
 12  2h2 sin  
KX
t=1
(vt;vt)
 12  2h2 sin  (v;v)
and noting that there exists  > 0 such that PKt=1(vt;vt)  (v;v), we also have
(E Bv;v) 
KX
t=1
 h2t sin t(vt;vt)(25)
  h2
KX
t=1
(vt;vt)
  h2 (v;v):
The inequalities (24) and (25) prove the theorem.  
7
Theorem 4. Let f2S T. Then hpC1jjfjjL2(D) L2( )  jjfjj2, where C1 is as in Theorem 3.
Proof. Using Theorem 1 we have jjfjj22 = ((E B)^f;(E B)^f). Now setting g = (E B) 12 ^f and
using Theorems 2 and 3 we have
C1h2  ((E B)g;g)(g;g) = ((E B)
^f;(E B)^f)
((E B)^f;^f) =
jjfjj22
jjfjj2L2(D) L2( )(26)
as required.  
Theorem 5. The inequality
C3h2  (Av;v)(v;v)  C4(27)
holds for all v2RMN, where C3 = C3(p) and C4 = C4(p).
Proof. We omit the proof. Su ce to say that it follows the proof for the analogous deterministic
problem as given in Elman et al. (2005).  
Note that Theorem 5 implies that the maximum eigenvalue of A is bounded above by a constant,
viz., C4.
2.6. Doubly Orthogonal Polynomials. We here establish two results concerning a basis for T
comprised of the so-called doubly orthogonal polynomials. The doubly orthogonal polynomials are
discussed in Babu ska et al. (2004). They are the set f 1;:::; Mg that satisfyZ
 
 k l =  kl;
Z
 
 r k l =  rk kl;(28)
where  kl is the Kronecker delta function.
An important consequence of the doubly orthogonal polynomials is that when they are used as a
basis for T the system matrix inx2.4 becomes block diagonal. This allows us to think of the discrete
problem as M decoupled (deterministic) steady-state di usion problems. As this is an important
point that will be used in the proof of the approximation property in x3.4 we state it explicitly in
the following theorem. This is then followed by a theorem concerning the source functions of these
decoupled problems that will also be used in proving the approximation property.
Theorem 6. If the  nite element approximation uhp is expanded using the doubly orthogonal poly-
nomials, viz.,
uhp =
MX
k=1
uk k;(29)
then each uk2S is the  nite element approximation to a steady-state di usion problem.
Proof. If f 1;:::; Mg is used for the computational basis then the system matrix in x2.4 will be
of the form
A =
2
64
A1
...
AM
3
75; [A
k]ij =
Z
D
bkr i r j;(30)
where
bk = c0 +
mX
r=1
p
 r rkcr(31)
where cr,  rk, r = 1;:::;m are from (2) and (28) respectively. This proves the theorem.  
8
Theorem 7. Let f 2S T and consider the expansion of the  nite element approximation uhp
is terms of the doubly orthogonal polynomials as given in (29). Then the source function fk of the
steady-state di usion problem to which uk is the solution satis es
jjfkjjL2(D)  jjfjjL2(D) L2( ):(32)
Proof. As f2S T we can expand it in the doubly orthogonal polynomials as
f =
MX
l=1
^fl l:(33)
Using f 1;:::; Mg as the computational basis and substituting (33) into (13) we have
[fk]i =
Z
 
Z
D
 MX
l=1
^fl l
 
 i k =
MX
l=1
Z
 
 l k
Z
D
^fl i =
Z
D
^fk i:(34)
This shows that fk = ^fk. Now
jjfjj2L2(D) L2( ) =
Z
 
Z
D
 MX
l=1
^fl l
 2
=
MX
l=1
Z
 
 2l
Z
D
^f2l(35)
=
MX
l=1
Z
D
^f2l =
MX
l=1
jjfljj2L2(D)  jjfkjj2L2(D);
for k = 1;:::;M, which proves the theorem.  
3. Multigrid
In this section we give a two-grid correction scheme for solving the system of equations given
in x2.4. This scheme varies the mesh parameter from grid to grid, i.e. there is a coarse grid and
a  ne grid, while the stochastic discretization parameter, p, is held constant. Thereby, the scheme
resembles that that would be applied to the regular deterministic problem. It is known that such a
scheme when applied to the deterministic problem will converge at a rate independent of the value
of the mesh parameter, h. We show that this is also the case for the stochastic problem. To show
this we follow a regular multigrid analysis, as given e.g. in Braess (2001) or Elman et al. (2005),
and show that a smoothing property and an approximation property hold. In order to establish the
approximation property we make use of the doubly orthogonal polynomials introduced inx2.6. Once
the convergence of the two-grid scheme has been shown to be independent of the mesh parameter it
follows, by an inductive argument, that the convergence of a multigrid algorithm, which applies the
two-grid algorithm recursively, is also independent of the mesh parameter.
3.1. Stationary Iteration. Central to the idea of multigrid is the understanding that certain
stationary iterations when applied to particular matrix problems tend to smooth the associated
error.
Consider a general matrix problem Au = f. Then the matrix splitting A = M N inspires the
stationary iteration
u(k+1) = M 1Nu(k) +M 1f(36)
= M 1(M A)u(k) +M 1f
= (I M 1A)u(k) +M 1f:
The matrix I M 1A is the iteration matrix of the method and in the context of multigrid is called
the smoother.
9
3.2. Two-grid Correction Scheme. Let T  L2( ) and S2h  Sh  H10 (D) be as de ned in
x2.3. Then de ning V 2h = S2h T and Vh = Sh T we have V 2h Vh H10 (D) L2( ). Finite
element formulations in Vh and V 2h give rise to matrix equations which we represent as Au = f
and  A u =  f respectively.
We now de ne a prolongation operator Ih2h: V 2h!Vh via natural inclusion, i.e. for v2h2V 2h,
Ih2hv2h = v2h. To see how Ih2h can be represented we note that any basis function  2hj of S2h can be
expanded in the basis functions of Sh, viz.,
 2hj =
NhX
i=1
pij hi; j = 1;:::;N2h:(37)
We de ne a matrix P using the coe cients above, i.e. [P]ij = pij. Now we have, for v2h2V 2h,
v2h =
N2hX
j=1
MX
k=1
v2hjk 2hj  k =
N2hX
j=1
MX
k=1
v2hjk
NhX
i=1
pij hi k(38)
=
NhX
i=1
MX
k=1
 N2hX
j=1
pijv2hjk
 
 hi k =
NhX
i=1
MX
k=1
[Pv2hk ]i hi k:
As v2h2Vh we also have the expansion
v2h =
NhX
i=1
MX
k=1
vhik hi k:(39)
Comparing (38) and (39) we see that [Pv2hk ]i = vhik or that Pv2hk = vhk. From this it follows that
if v2h is the coe cient vector of v2h in V 2h, then (I P)v2h is the coe cient vector of v2h in Vh.
(Here I is an M M identity matrix.) We therefore call I P the prolongation matrix and introduce
the notation P = I P.
We next de ne a restriction operator I2hh : Vh ! V 2h such that the corresponding restriction
matrix R satis es R=PT (or equivalently R= I R where R = PT). That is to say, that if I2hh
maps vh2V 2h to v2h2V 2h and vh and v2h are the respective coe cient vectors of these functions,
then v2h = Rvh = PTvh. With the prolongation and restriction operators related in this way we
have the desirable relationships  f =Rf and  A =RAP.
Using these de nitions we have the following algorithm for a two-grid iterative correction scheme.
choose initial guess u
for i = 0;1;:::
for j = 1 : k
u (I M 1A)u +M 1f
end
 r =R(f  Au)
solve  A e =  r
u u +P e
end
The success of this algorithm necessarily depends on how well the smoother works and how well the
functions are passed between the coarse and  ne grids.
3.3. Convergence of Two-Grid Correction Scheme. We wish to establish that the two-grid
algorithm given in x3.2 converges and that the contraction rate is independent of h. This can be
10
shown to be true providing the smoothing property and the approximation property are satis ed, as
is shown in the following theorem.
Theorem 8. Providing the smoothing property,
jjA(I M 1A)kyjj2   (k)jjyjjA 8y2RMNh;(40)
with n(k)!0 as k!1, and the approximation property,
jj(A 1 P  A 1R)yjjA C5jjyjj2 8y2RMNh;(41)
with C5 = C5(p), are satis ed, then, providing k is su ciently large, the two-grid algorithm given in
x3.2 converges and the contraction rate is independent of h.
Proof. It can be shown that the error associated with the two-grid algorithm obeys the recursive
relationship
e(i+1) = (A 1 P  A 1R)A(I M 1A)ke(i):(42)
Hence,
jje(i+1)jjA =jj(A 1 P  A 1R)A(I M 1A)ke(i)jjA(43)
 C5jjA(I M 1A)ke(i)jj2
 C5 (k)jje(i)jjA:
Since  (k)!0 as k!1there exists some minimal number of smoothing steps such that C5 (k) <
1.  
The proof that the smoothing property holds is dependent on the choice of smoother, i.e. the
value of M. For the case of M =  I,  2R, which choice corresponds to Richardson?s iterative
method, the proof follows that given in Braess (2001) and Elman et al. (2005). We will prove that
the approximation property holds in x3.4.
3.4. Approximation Property. We here wish to show that the approximation property given in
(41) is satis ed.
Theorem 9. For the problem under consideration, the approximation property given in Theorem 8
holds.
Proof. Given y2RMNh we can  nd some f2Sh T such that y = f. Let uhp and u2h;p be the  ne
and coarse grid solutions respectively with coe cient vectors u = A 1f and  u =  A 1 f =  A 1Rf.
Then we have
jj(A 1 P  A 1R)yjj2A =jju P ujj2A = (u P u;u P u)A(44)
= a(uhp Ih2hu2h;p;uhp Ih2hu2h;p)
= a(uhp u2h;p;uhp u2h;p):
Utilizing the fact that c is bounded above P-a.e. by  and expanding uhp and u2h;p using the doubly
orthogonal polynomials we have
jj(A 1 P  A 1R)yjj2A  
Z
 
Z
D
  
  r
 MX
k=1
uhk k 
MX
k=1
u2hk  k
   
  
2
:(45)
11
Now using the fact that the doubly orthogonal polynomials are orthonormal we get
jj(A 1 P  A 1R)yjj2A  
MX
k=1
Z
 
 2k
Z
D
jr(uhk u2hk )j2(46)
=  
MX
k=1
Z
D
jr(uhk u2hk )j2
  
MX
k=1
jjuhk u2hk jj2H1(D):
From Theorem 6 we know that uhk2Sh, u2hk 2S2h, k = 1;:::;M, are  nite element approximations
to (deterministic) steady-state di usion problems. Let uk 2 H10 (D), k = 1;:::;M, be the weak
solutions to these problems. Furthermore, let fk2Sh, k = 1;:::;M, be the source functions of the
associated problems. Then there exists constant Cd > 0 such that jju uhkjjH1(D)  CdhjjfkjjL2(D)
for all k. Therefore,
jj(A 1 P  A 1R)yjj2A  
MX
k=1
(jjuk uhkjjH1(D) +jjuk u2hk jjH1(D))2(47)
  
MX
k=1
(CdhjjfkjjL2(D) + 2CdhjjfkjjL2(D))2
=  
MX
k=1
(3CdhjjfkjjL2( ))2:
Now using Theorems 4 and 7 we have
jj(A 1 P  A 1R)yjj2A  
MX
k=1
(3CdhjjfjjL2(D) L2( ))2(48)
  M
 3C
dp
C1jjfjj2
 2
=  M
 3C
dp
C1jjyjj2
 2
which completes the proof.  
3.5. Extension to Multigrid. The two-grid correction scheme given in x3.2 only contains pre-
smoothing. In practice post-smoothing is often also applied. In the numerical experiments given in
x4 post-smoothing is applied. We have neglected post-smoothing in the preceding analytic argument
in order to keep things a little simpler. It can be shown, though we omit the details here, that the
two-grid correction scheme with post-smoothing also converges with a contraction factor independent
of the spatial mesh parameter.
Recursively applying the two-grid correction scheme, gives rise to a multigrid scheme. A number
of variations are possible, see, for example, Briggs, Henson & McCormick (2000). That multigrid
converges with a contraction factor independent of the spatial mesh parameter can be established
by an inductive argument once the two-grid scheme has been shown to converge with a contraction
factor independent of the spatial mesh parameter. This inductive argument will be no di erent for
the stochastic problem than for the analogous deterministic problem and is discussed e.g. in Braess
(2001) and Elman et al. (2005).
12
4. Numerical Experiments
We now perform some numerical experiments to provide practical support for the theoretical
results obtained in x3. The model that we use follows that given in Ghanem & Spanos (1991)
and Xiu & Karniadakis (2002). We would also like to direct the reader?s attention to Ma^itre,
Knio, Debusschere, Najim & Ghanem (2003) where multigrid was applied to stochastic steady
and unsteady di usion problems and in which the conclusions reached are in agreement with the
conclusions that we reach in x4.6 and x4.7.
4.1. Model Problem. We take the spatial domain to be D = ( 1;1)2 and consider the determin-
istic source function f = 1.
To construct the di usion coe cient we consider a process with mean function c0(x), constant
variance  , and covariance function r(x;y). Such a process will have a Karhunen-Lo eve expansion
of the form
c(x;!) = c0(x) +
1X
k=1
p
 kck(x) k(!)(49)
where ( k) is a sequence of uncorrelated and identically distributed random variables, and ( k) and
(ck) can be computed by solving the eigenvalue equation
Z
D
r(x;y)ck(x) dx =   kck(y):(50)
If need be we make the further assumption that ( k) is a sequence of independent random variables.
The sequence ( k) is ordered so as to be non-increasing.
For computational purposes we need a  nite term expansion so we approximate (49) by
c(x;!) = c0(x) +
mX
k=1
p
 kck(x) k(!)(51)
where  k, ck, k = 1;:::;m, still satisfy (50). From the modeling perspective the replacement
of the in nite expansion with the  nite expansion is justi ed providing ( k) decays rapidly. We
demonstrate the decay of ( k) in x4.2 where we discuss the covariance function that we will use.
We will consider two cases for the distributions of the random variables  k, k = 1;:::;m. Inx4.6
we take  k, k = 1;:::;m, to be uniformly distributed on ( 1;1) with c0(x) = 10 and  = 1=3. In
x4.7 we take  k, k = 1;:::;m, to be normally distributed with c0(x) = 10 and  = 0:01.
For information on the Karhunen-Lo eve expansion the reader is referred to Lo eve (1994). And for
a discussion of its use in mathematical modeling the reader is referred to Ghanem & Spanos (1991).
4.2. Exponential Covariance. We consider for the covariance function of the di usion coe cient
the exponential covariance function given by
r(x;y) =  e 1bjx1 y1j 1bjx2 y2j(52)
where x = (x1;x2), y = (y1;y2)2D. The constant b is called the correlation length and will a ect
the decay of ( k), a larger value producing faster decay.
Analytic expressions for ( k) and (ck) can be obtained for this choice of covariance function.
The derivation is given in Ghanem & Spanos (1991). We here will simply state the results. Given
D = ( a;a)2, let
 1k = 2=b 2
k + (1=b)2
; k = 1;2;:::(53)
13
and
c1k(x) = cos( kax)q
a+ sin(2 ka)2 
k
; k = 1;3;:::;(54)
c1k(x) = sin( kax)q
a sin(2 ka)2 
k
; k = 2;4:::;(55)
where c1k: ( a;a)!Rand ( k) > 0 is an increasing sequence satisfying the transcendental equations
(1=b)  k tan( ka) = 0;  k + (1=b) tan( ka) = 0:(56)
Then
 k =  1i 1j; ck(x) = c1i(x1)c1j(x2);(57)
where the ordering is such that ( k) forms a non-increasing sequence. We note that ( 1k) and (c1k)
are in fact the eigenvalues and eigenfunctions of the equivalent 1-dimensional problem.
We give these analytic expressions of ( k) and (ck) so that the reader may better see the nature
of the di usion coe cient. In practice, however, the covariance function may be such that analytic
expressions for ( k) and (ck) are not available and (50) will need to be solved numerically. In fact,
in the numerical experiments given inx4.6 andx4.7 we solved (50) numerically. This is not a trivial
task as it involves solving a generalized eigenvalue problem in which the right hand side matrix, in
this context sometimes called the Galerkin covariance matrix, is in general dense. This can be done
e ciently by using the Lanczos algorithm together with fast methods for matrix-vector products,
as discussed in Eiermann, Ernst & Ullmann (2005). This issue is also discussed in Karniadakis, Su,
Xiu, Lucor, Schwab & Todor (2005).
For further discussion on appropriate choices of covariance functions see Ghanem & Spanos (1991)
and Xiu & Karniadakis (2002).
4.3. Matrix Expansions. With c and f de ned as inx4.1 we  nd that the matrices A and f given
in x2.4 have the form
A = G0 A0 +
mX
k=1
p
 kGk Ak;(58)
f = g0 f0(59)
where, de ning  0 = 1,
[Gk]ij =
Z
 
 k i j; [g0]i =
Z
 
 i;(60)
[Ak]ij =
Z
D
ckr i r j; [f0]i =
Z
D
 i:(61)
We note that the matrices Ak will be sparse due to the choice of spatial basis functions for S. These
will form the blocks of the system matrix, A. The block structure will be determined by the matrices
Gk and so it is important that the basis functions for T are chosen so that these too are sparse.
4.4. Polynomial Chaos. There are a number of choices on how to construct T. As mentioned in
x2.3 we follow the methodology given in Xiu & Karniadakis (2002) which generalizes the method
of polynomial chaos as given in Ghanem & Spanos (1991). The functions thus chosen to be the
computational basis for T are those from the Askey scheme of hypergeometric polynomials that are
orthogonal to with respect to the probability measure, i.e.Z
 
 k l = dk kl;(62)
14
0 2 4 6 8 10 12 14 16
0
2
4
6
8
10
12
14
16 nz = 55
m = 4, p = 2
0 2 4 6 8 10 12 14 16
0
2
4
6
8
10
12
14
16 nz = 55
m = 2, p = 4
Figure 1. Block structure of system matrix A.
and are of degree p or less. It follows directly that the matrix G0 will be diagonal. Furthermore,
the matrices Gk, k = 1;:::;m, will turn out to be sparse due to the three-term recurrence relation
that all orthogonal polynomials satisfy. The resultant block structure of A is shown in Figure 1 for
two choices of m and p. We emphasize that the solid block shown in the  gure represents a sparse
matrix having the struture of a sti ness matrix arising from a deterministic steady-state di usion
problem, as discussed in x4.3.
4.5. Spatial Mesh. As mentioned in x4.1 the spatial domain is taken to be D = ( 1;1)2. For the
triangulation of D we will use uniform meshes consisting of an underlying grid of n n squares each
of which is further subdivided into two equal triangles.
4.6. Multigrid for Di usion with Uniform Distributions. We now let  k, k = 1;:::;m, be
uniformly distributed on  = ( 1;1). Therefore,  = ( 1;1)m and dP = d!=2m. We also set
 = 1=3.
Given the nature of  k, ck, k = 1;:::;m, as given inx4.2, and the form of the di usion coe cient
as given in x4.1, it can be shown that (6) will be satis ed providing the mean is su ciently large.
We here take c0(x) = 10.
Applying the generalized polynomial chaos method as described inx4.4 the basis of T will be the
set of m-dimensional Legendre polynomials of degree p or less.
Now multigrid is applied. A full V-cycle is used with an n n  nest mesh and a 2 2 coarsest
mesh. For the smoother we use the damped Jacobi method with the damping parameter set to 2/3.
Three pre-smoothing and three post-smoothing iterates are carried out. The iterations stop when
the relative residual reaches a tolerance of 10 6. Table 1 shows the number of iterations required for
convergence for varying values of m, n, and p. The results clearly support the theoretical conclusion
that the contraction rate of the multigrid algorithm is mesh independent. The table also indicates
that the method is apparently insensitive to the parameters m and p.
4.7. Multigrid for Di usion with Normal Distributions. We now let  k, k = 1;:::;m,
be normally distributed with zero mean and variance  . Now we have  = Rm and dP =
e !2=(2 )=(2  )m=2. We take c0(x) = 1.
Note that the di usion coe cient as de ned inx4.1 will now fail to satisfy condition (6) no matter
what the choice of  . However, we have reason to believe that the theory still applies. We give here
only a heuristic argument. Given a su ciently small variance the probability of c being outside of
two positive bounds becomes negligibly small. That is to say, that if the normal distributions were
replaced by similar distributions that looked like the normal distributions with their tails cut o 
15
n = 4 m = 1 m = 2 m = 3 m = 4
p = 1 6 6 6 6
p = 2 6 6 6 6
p = 3 6 6 6 6
p = 4 6 6 6 6
n = 8 m = 1 m = 2 m = 3 m = 4
p = 1 7 7 7 7
p = 2 7 7 7 7
p = 3 7 7 7 7
p = 4 7 7 7 7
n = 16 m = 1 m = 2 m = 3 m = 4
p = 1 7 7 7 7
p = 2 7 7 7 7
p = 3 7 7 7 7
p = 4 7 7 7 7
n = 32 m = 1 m = 2 m = 3 m = 4
p = 1 7 7 7 7
p = 2 7 7 7 7
p = 3 7 7 7 7
p = 4 7 7 7 7
Table 1. Number of iterations required for multigrid to converge for di usion
de ned via uniform distributions.
so as to ensure that c satis es (6), then the di erence would not be noticed computationally. We
emphasize that we have not pursued this reasoning analytically. We have found that su ciently
small variance results in positive de nite systems that yield sensible results. We note that this
problem has been tackled in Ghanem & Spanos (1991) and Xiu & Karniadakis (2002). We take
 = 0:01.
Applying the generalized polynomial chaos method as described inx4.4 the basis of T will be the
set of m-dimensional generalized Hermite polynomials of degree p or less.
Now multigrid is applied. A full V-cycle is used with an n n  nest mesh and a 2 2 coarsest
mesh. For the smoother we use the damped Jacobi method with the damping parameter set to 2/3.
Three pre-smoothing and three post-smoothing iterates are carried out. The iterations stop when
the relative residual reaches a tolerance of 10 6. Table 2 shows the number of iterations required for
convergence for varying values of m, n, and p. The results support the theoretical conclusion that
the contraction rate of the multigrid algorithm is mesh independent. The table also indicates that
the method is apparently insensitive to m and only slightly sensitive to p.
References
Babu ska, I., Tempone, R. & Zouraris, G. E. (2004), ?Galerkin  nite element approximation of stochastic elliptic partial
di erential equations?, Siam J. Numer. Anal. 42, 800{825.
Braess, D. (2001), Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics, second edn, Cam-
bridge University Press.
Briggs, W. L., Henson, V. E. & McCormick, S. F. (2000), A Multigrid Tutorial, second edn, SIAM.
Eiermann, M., Ernst, O. G. & Ullmann, E. (2005), Computational aspects of the stochastic  nite element method, in
?Proceedings of Algoritmy?.
16
n = 4 m = 1 m = 2 m = 3 m = 4
p = 1 6 6 6 6
p = 2 7 7 7 7
p = 3 7 7 7 7
p = 4 7 7 7 7
n = 8 m = 1 m = 2 m = 3 m = 4
p = 1 8 8 8 8
p = 2 8 8 8 8
p = 3 9 9 9 9
p = 4 10 10 10 10
n = 16 m = 1 m = 2 m = 3 m = 4
p = 1 8 8 8 8
p = 2 8 8 8 8
p = 3 9 9 9 9
p = 4 9 10 10 10
n = 32 m = 1 m = 2 m = 3 m = 4
p = 1 7 7 8 8
p = 2 8 8 8 8
p = 3 8 8 9 9
p = 4 9 9 9 9
Table 2. Number of iterations required for multigrid to converge for di usion
de ned via uniform distributions.
Elman, H. C., Silvester, D. J. & Wathen, A. J. (2005), Finite Elements and Fast Iterative Solvers: with Applications
in Incompressible Fluid Dynamics, Oxford University Press.
Ghanem, R. G. & Spanos, P. D. (1991), Stochastic Finite Elements: A Spectral Approach, Springer-Verlag.
Karniadakis, G., Su, C.-H., Xiu, D., Lucor, D., Schwab, C. & Todor, R. (2005), Generalized polynomial chaos solution
for di erential equations with random inputs, Technical Report 2005-01, Eidgen ossische Technische Hochschule
Z urich.
Lo eve, M. (1994), Probability Theory II, Springer.
Ma^itre, O. L., Knio, O., Debusschere, B., Najim, H. & Ghanem, R. (2003), ?A multigrid solver for two-dimensional
stochastic di usion equations?, Comput. Methods Appl. Mech. Engrg. 192, 4723{4744.
Treves, F. (1967), Topological Vector Spaces, Distributions and Kernels, Academic Press.
Xiu, D. & Karniadakis, G. E. (2002), ?Modeling uncertainty in steady state di usion problems via generalized poly-
nomial chaos?, Comput. Methods Appl. Mech. Engrg. 191, 4927{4948.