COMPUTATIONAL METHODS IN APPLIED MATHEMATICS, Vol.1(2001), No.3, pp.

222–241
c Institute of Mathematics of the National Academy of Sciences of Belarus
°

ON THE REPRESENTATION OF FUNCTIONS

AND FINITE DIFFERENCE OPERATORS ON
ADAPTIVE DYADIC GRIDS
P.W. HEMKER
CWI
P.O. Box 94079, 1090 GB Amsterdam, The Netherlands
E-mail: [email protected]

F. SPRENGEL
Institute for Algorithms and Scientific Computing
GMD German National Research Center for Information Technology
Schloss Birlinghoven, D–53754 Sankt Augustin, Germany
E-mail: [email protected]

Abstract — In this paper we describe methods to approximate functions and dif-

ferential operators on adaptive sparse (dyadic) grids. We distinguish between several
representations of a function on the sparse grid and we describe how finite difference
(FD) operators can be applied to these representations.
For general variable coefficient equations on sparse grids, genuine finite element (FE)
discretizations are not feasible and FD operators allow an easier operator evaluation
than the adapted FE operators. However, the structure of the FD operators is complex.
With the aim to construct an efficient multigrid procedure, we analyze the structure of
the discrete Laplacian in its hierarchical representation and show the relation between
the full and the sparse grid case. The rather complex relations, that are expressed by
scaling matrices for each separate coordinate direction, make us doubt about the pos-
sibility of constructing efficient preconditioners that show spectral equivalence. Hence,
we question the possibility of constructing a natural multigrid algorithm with optimal
O(N ) efficiency.
We conjecture that for the efficient solution of a general class of adaptive grid problems
it is better to accept an additional condition for the dyadic grids (condition L) and to
apply adaptive hp-discretization.
2000 Mathematics Subject Classification: 65D15; 65N06; 65N30; 65N55.

Keywords: sparse grid, dyadic grid, adaptive methods.

1. Introduction
When applied to d-dimensional problems (d > 2), all numerical methods using regular
rectangular grids have one problem in common: the curse of dimension. This means that
if one refines the grid — for instance, by repeatedly halving the mesh size — the number
of grid points grows exponentially with the dimension, i.e., like N d in Rd where N = O(2n )
 Finite difference operators on dyadic grids 223

denotes the number of points per direction. One way out seems to be the use of sparse grids,
where the number of points only grows like N (log N )d−1 . Under certain conditions on the
mixed derivatives of the approximated function, for sparse grids the approximation accuracy
degrades only by a logarithmic factor compared with the accuracy achieved on the regular
grid with the same N , see [17, 20, 42].
If we study the efficient representation of functions on a grid in more dimensions, i.e., if
we relate the accuracy of the approximation to the number of degrees of freedom (number
of data points) used, then it becomes clear that the assumed regularity of the approximated
function is the essential issue. For sufficiently smooth functions, more efficient methods
are obtained by higher order approximation. In this light, sparse grid approximation can
be considered as a method where a particular regularity assumption (viz., bounded mixed
derivatives) is exploited, that is intermediate between usual assumptions of C k -regularity.
More efficient than the usual sparse grid approximation, by piecewise (multi-) linear func-
tions, is approximation by a higher order. The same argument holds for higher order sparse
grid approximations.
More efficient representations are further obtained by applying mesh refinement not over
the whole domain, but only in those regions where a function has no sufficient regularity
to apply the higher order approximation. So we arrive at the use of adaptive grids. By
simultaneous use of many regular, overlapping grids, for which each time the mesh width is
halved in one of the coordinate directions, and by using only those parts of the grids that
are needed for an efficient representation, we arrive at the adaptive dyadic grid structure.
This structure may possibly yield an optimally efficient algorithm for the solution of partial
differential equations. In fact, both the usual sparse and the regular full grids are special
cases for this structure.
Sparse grids and related methods already have a long tradition in numerical quadrature
and approximation theory (see, e.g., [7, 22, 34, 39]). During the last decade, since sparse grid
methods have been introduced for the numerical treatment of elliptic boundary problems by
Zenger [42], several authors have used sparse grids for solving partial differential equations
with finite elements [4–6, 12, 17, 26, 29, 31], finite volumes [18] or spectral methods [23].
Beside this, sparse grids are nowadays used for high-dimensional approximation [3, 37] and
numerical integration of functions needed, e.g., in financial mathematics [10,28,41], in wavelet
theory [8, 35] or in visualization [38].
Since the efficient evaluation of genuine variable coefficient 3-dimensional finite element
operators is impossible, modified Galerkin-type discretizations have been proposed [1, 9,
30]. However, generally these discretizations loose the typical Galerkin nesting property
and, e.g., symmetric differential operators turn into asymmetric difference operators. Thus,
the proposed Galerkin methods are essentially not more attractive than a finite difference
method. Hence, recently also finite difference methods for sparse grids have been developed
[13, 16, 32, 33]. Note that consistent finite difference operators on sparse grids cannot be
chosen as simple stencils involving the nearest neighbors of a point, cf. [32; Appendix].
In [13, 32], the whole machinery for finite differences on sparse grids (which needs linear
combination, multiplication, and approximate differentiation) is provided. In Section 3, we
show a simple way of presenting the finite difference discretization, and we show a direct
construction of the second derivative from the hierarchical coefficients, by recursion.
In Section 4, we restrict ourselves to the special case of Poisson’s equation and regular
sparse grids in order to analyze the finite difference operator in detail. The resulting matrix
is ill-conditioned. So, it takes (even in the Jacobi-preconditioned version [21, 32]) many
224 P.W. Hemker and F. Sprengel

iterations for an iterative solver (like BiCGStab) to obtain a solution. Our purpose is to
construct Galerkin relations which may motivate efficient multilevel algorithms. An attempt
to obtain such an algorithm is demonstrated. However, like [32], we did not succeed in finding
a method with a level-independent convergence rate. The Galerkin relations appear to be
complex and we doubt if spectral equivalence can be realized. Hence, finally, we conclude
that it may be practical to restrict the adaptive structure further in order to obtain a dyadic
hp-adaptive structure satisfying ‘condition L’. For such a structure, in which a unique ‘finest’
discretization is available, Galerkin-type methods give no special problems.

2. Representation of ASG functions

2.1. Basic notation
To be able to describe adaptive dyadic or sparse grid (ASG) function representations, we
first summarize some necessary notation. For background see, e.g., [20].
• Domain Ω ⊂ Rd , with coordinates xj , j = 1, . . . , d.
• Multi-integer: m = (m1 , m2 , . . . , md ) ∈ Nd0 ,

– 0 = (0, 0, . . . , 0),
– e = (1, 1, . . . , 1),
– ej = (. . . , 0, 1, 0, . . .), the j-th unit vector,
P
– |m| = dj=1 mj ,
Q
– |||m||| = dj=1 mj ,
– m < n ⇔ mj < nj ∀j = 1, 2, . . . , d.

• Dyadic mesh Ωk , k > 0, mesh with dyadic mesh-width hk ,

– mesh-width: hk = (hk1 , hk2 , . . . , hkd ) with hki = 2−ki ,

– dyadic grid: Ω+ k = {xk,j | xk,j = jhk = (j1 hk1 , j2 hk2 , . . . , jd hkd )} ∩ Ω,
i.e., the grid of nodal
S points in the mesh Ωk ,
+ +
– sparse grid: Ω` = |k|=` Ωk ,
– a point with |||m||| odd is called a hierarchical point,
– with “|||j||| odd” we mean: for all i = 1, 2, . . . , d, either ji is an odd integer, or ki = 0.
Q mj
• Derivatives: Dm = dj=1 ∂∂xj .
• Univariate hat function: ϕ(x) = max(0, 1 − |x|).
• Univariate Haar function: η(x) = 1 for 0 < x < 1, and η(x) = 0 for x < 0 or x > 1.
Q
• Basis hat functions: ϕk,j (x) = di=1 ϕ(xi /hki − ji ).
Q
• Basis Haar functions: ηk,j (x) = di=1 η(xi /hki − ji ).
• Space of piecewise d-linear functions on Ωk : Vk = span{ϕk,j | j ∈ Zd , xk,j ∈ Ω+
k }.
• Space of hierarchical surpluses on Ωk :
Wk = span{ϕk,j | |||j||| odd, j ∈ Zd , xk,j ∈ Ω+
k}.

Without loss of generality, we assume here that k = 0 yields “the coarsest grid”.
 Finite difference operators on dyadic grids 225

2.2. The hierarchical (H-) and lattice (L-) condition

Given a continuous function u ∈ C(Ω), we can approximate it by the interpolant un ∈ Vn
by means of interpolation on the grid Ω+
n , i.e.,

un (xn,j ) = u(xn,j ) ∀xn,j ∈ Ω+

n .

Obviously, the function un on Ωn is given by

X
un = an,j ϕn,j , (1)
j

where an,j = u(xn,j ). The error of approximation is wellknown. However, in contrast to

the classical approximation we are not interested in approximation for a fixed n, but in the
approximation on (the union of) a number of grids Ω+ n.
We can make the approximation (1) for all grids Ω+ n with n > 0. It can be made
arbitrarily accurate for n large enough, but the number of degrees of freedom increases
geometrically with |n|. Therefore, in practice we select a ‘smallest’ n such that an accuracy
criterion is satisfied. Notice that keeping the representations in all coarser Vk (all Vk , 0 6
k 6 n) does not take essentially more coefficients than the representation on the finest grid
in Vn alone.
In order to obtain an efficient approximation, it may be useful to distinguish among
different subregions where we make the finest approximation of u in a different Vn . We
make full and efficient use of the system {Vn | n ∈ Nd0 }, by in principle approximating a
given function u ∈ C(Ω) in all {Vn | n ∈ Nd0 }, but using in practice only those coefficients
that contribute to a sufficiently accurate representation. This implies that, in practice, the
function u can be represented in a particular Vn only on part of the domain Ω. To introduce
structure in the family of approximating basis functions {ϕn,j }, we introduce the following
conditions H and L, which respectively introduce locally a hierarchical and a lattice structure.
Condition H: If a basis function ϕn,j is used in the representation (1), all corresponding
coarser basis functions, i.e., functions ϕk,i for which supp(ϕk,i ) ⊃ supp(ϕn,j ), are also used
for the approximation.
Condition L: If two basis functions ϕn1 ,j1 and ϕn2 ,j2 with supp(ϕn1 ,j1 )∩supp(ϕn2 ,j2 ) 6= ∅
are used in the representation (1), then the corresponding finer basis function, i.e., a function
ϕn3 ,j3 for which supp(ϕn3 ,j3 ) = supp(ϕn1 ,j1 )∩supp(ϕn2 ,j2 ), is also used for the approximation.
Condition H is a minimal structure required for sparse grids. Much more structure
appears in the partially ordered family of grids {Ωn } used for the approximation if an
approximation satisfies additionally the condition L. Then the grids in the (local) approx-
imation form a lattice. Condition L holds for the grids used in multigrid methods with
semi-coarsening [24, 25]. It guarantees that all coarse grid approximations can be obtained
by restriction from a unique finest grid. Typically, sparse grids do not satisfy condition L.

2.3. E-, C-, D- and H-representation

We call the representation of the approximation of a function u ∈ C(Ω) by a collection
of such (partial) approximations in the family of spaces {Vn } the nodal representation or
the E-representation of the approximation. This E-representation requires the coefficients
an,j = u(xn,j ), corresponding to grid points xn,j , to be equal on the different grids Ω+n at
coinciding grid points xn,j . Thus, because points from coarser grids coincide with those from
226 P.W. Hemker and F. Sprengel

finer ones, a certain consistency is required (and a redundancy exists) in the E-representation
of an approximation.
During the computation in an approximation process, the representations of the approx-
imations on all different grids Ω+k do not necessarily always satisfy the consistency condition
required for the E-representation. In that case an approximation exists of the form of (1)
on each separate grid Ω+ n , and the approximation on the whole system is not uniquely de-
termined. Such a representation is called the collateral or C-representation of a (nonunique)
approximation. There, for different n, the approximations un (x) do not necessarily coincide
at corresponding grid points xn,j .
Another way of representing approximations on the family of grids {Ω+ n } is to partition
the approximation over the different grids. Then, instead of (1) the approximation reads
XX
uh = an,j ϕn,j . (2)
n j

In this case, of course, the set of coefficients {an,j } always determines a unique function uh .
An approximation in this form is called a distributed or D-representation. However, for a
given function uh , now the coefficients {an,j } are not uniquely determined because the {ϕn,j }
are linearly dependent.
One way to select a special unique D-representation is by choosing the coefficients an,j
such that an,j 6= 0 only for those (n, j) for which |||j||| is odd. This implies that an,j = 0 except
for a pair (n, j) for which Ω+n is the coarsest grid which contains the nodal point xn,j . This
representation X
uh = an,j ϕn,j (3)
(n,j),|||j||| odd

is called the H-representation because it represents the approximation in the hierarchical

basis © ª
ϕn,j | n ∈ Nd0 , j ∈ Zd , |||j||| odd, xn,j ∈ Ω+
n , (4)
and the part of uh in
Wn = span{ϕn,j | j ∈ Zd, |||j||| odd, xn,j ∈ Ω+
n}

is the hierarchical contribution from the grid Ω+

n to the approximation. We notice that
d
X X
Vn = Wn + Vn−ej = Wm ,
j=1 6 6
0 m n

and that the sparse grid space is defined by

X
VL = Wm .
6 6
0 |m| L

When u is interpolated at the nodal points xn,j , the hierarchical coefficients an,j in
X
u(xn,j ) = an,j ϕn,j (xn,j ) (5)
(n,j),|||j||| odd

are determined by (cf. [20])

Yd · ¸
1 1
an,j = − , 1, − u(jhn ) , (6)
2 2 hn
i=1 i ei
 Finite difference operators on dyadic grids 227

£ ¤
where − 21 , 1, − 12 hn ei denotes the difference stencil for the mesh-size hni in the i-th coordi-
i
nate direction. Notice that this expression is well defined for each odd j because Condition
H requires that all hi -neighbors are nodal points in the approximation. Another expression
for the coefficient an,j is found in the following lemma, see [20].

Lemma 1. Let u∈C e+m , for a given m with 0 6 m 6 e. Then, for each ϕn,j ∈Wn , we
have
Yd · ¸
1 1
an,j = − , 1, − u(jhn )
i=1
2 2 hn ei
Zi
(7)
|e+m| −d−|n| e+m e−m n
= (−1) 2 D u(x) D ϕ(2 x − j)dΩ.
Ω

2.4. Transformation of representations and data structure used

From the above, it is clear that each H-representation is a D-representation and each E-
representation a C-representation. For piecewise d-linear functions, it is often described
[11–13] how a pyramid algorithm can be used to convert an E-representation to a H-
representation, and vice versa. Such a conversion can be executed in O(N ) operations,
where N is the total number of coefficients (degrees of freedom). The transformation from
a D-representation to an H-representation is equally straightforward.
The E-, H-, D-, and C-representations can also be used for piecewise constant functions,
and — because of the tensor product structure — discrete function representations can be
combined in the different coordinate directions. For example , a discrete function can be
piecewise constant in one and piecewise linear in the other coordinate directions. Also, for
the piecewise constant functions, efficient pyramid conversion algorithms exist between the
different (H-, D-, E-) representation styles. In this case, one should opt for either left- or
right-continuity at the discontinuities in the representation.
The data structure used to implement all the above possibilities of an adaptive (sparse)
grid representation can be efficient and relatively simple. For the d-dimensional case (d =
1, 2, 3), we use the data structure BASIS3 [19,27] that takes the ‘patch’ Pn,j as an elementary
entity. This Pn,j takes all information related to a right-open left-closed cell:
3
Y £ ¢
jk 2−nk , (jk + 1)2−nk .
k=1

This implies that there exist as many patches in the data structure as there are points used
in the description of the approximation. The patches are related to each other by means
of pointers in an intertwined tree structure, where each patch has at most 15 pointers to
related patches (3 fathers, 6 neighbors, and 6 kids). The data structure is symmetric with
respect to any of the coordinate directions.

3. Evaluation of difference operators to ASG functions

Although finite element (FE) discretization of a PDE on a sparse grid is feasible for a constant
coefficient problem in two dimensions, finite elements for higher-dimensional problems and
variable coefficients give rise to problems. The difficulty arises because, with the hierarchical
228 P.W. Hemker and F. Sprengel

basis (4) for test and trial space, the computational complexity of the evaluation of the
discrete operator becomes too large. This is caused by the fact that the intersection of the
supports of an arbitrary trial and test function is much smaller than the supports of these
functions themselves. This has as a consequence that the advantage of sparse grids is lost if
the FE discrete operator is evaluated.
The alternative is to essentially modify the FE discretization as in [1, 9, 30] or, as it was
already proposed in [13, 32], to use a finite difference discretization. To this end we should
be able to apply (approximate) differentiation to discrete representations of approximations
as described in Section 2. Then, the application of linear difference operators approximating
the linear differential operator
X ∂ µ ∂
¶ X
∂
Aij (x) + Bi (x) + C(x) (8)
i,j
∂xi ∂xj i
∂xi

comes down to the construction of linear combinations of ASG functions, and pointwise mul-
tiplication and approximate differentiation of such functions. In any of the representations
from Section 2.3 the construction of linear combinations is directly computed by linear com-
bination of the coefficients. Pointwise multiplication is only possible in the E-representation,
in which the function values at grid points are directly available. Below we describe differ-
entiation, which requires some more attention, distinguishing between the evaluation of first
and second order derivatives. Essentially, FD discretizations are also found in [13, 32] but
we take a slightly different point of view which simplifies the description.

3.1. First order derivatives

For a piecewise d-linear ASG function, the derivative ∂x∂ i uh (x) is well defined almost every-
where. Written in D-representation (2), its derivative is simply described by
∂ X
uh (x) ≡ Di uh (x) = an,j Di ϕn,j (x)
∂xi n,j

X d
Y
= an,j Di ϕ(2nk xk − jk )
n,j k=1

X d
Y
ni ni
= an,j (η(2 xi − ji + 1) − η(2 xi − ji )) ϕ(2nk xk − jk ) .
n,j k=1,k6=i

This, again, is a function in D-representation, piecewise constant in the i-direction and

piecewise linear in the other directions. It can be described by coefficients associated with
nodal points if we decide to choose either a left- or a right-continuous representation.

3.2. Second order derivatives

The computation of second order derivatives of piecewise d-linear ASG functions, Di2 uh (x),
seems to be less obvious because second derivatives of the piecewise d-linear functions vanish
almost everywhere on Ω. Nevertheless the approximation of the second order derivatives is
useful and can easily be derived from the representation of uh .
In order to approximate Di2 uh , we first construct the representation which is the H-repre-
sentation in the i-th coordinate direction and the E-representation in the other coordinate
 Finite difference operators on dyadic grids 229

directions, i.e., we apply the pyramid algorithm only in the i-th direction. This implies,
cf. (7),
· ¸
1 1
an,j = − , 1, − u(jhn )
2 2 hn ei
Z i (9)
−1−ni
= −2 Di2 u(x) ϕ(2ni xi − ji )dxi ,

where xk = 2−nk jk for k 6= i. The difference from [13, 32] is that we use the hierarchi-
cal coefficients in the i-th direction directly instead of forming second differences of an
E-representation. This gives us the possibility of describing the approximation of the second
derivatives in hierarchical representation explicitly in Section 4.
For the representation of the second derivative, we see that the coefficients in expression
(2) are given by (7). It follows that the hierarchical coefficient an,j with respect to the i-th
coordinate direction corresponds with a measure of Dx2i uh (x) in the hn -neighborhood of xn,j .
Considering the i-th coordinate in (9), writing h = 2−n , and omitting higher order terms in
h, we see
Z
h
an,j =− Di2 u(x)ϕ(x/h − j) dx
2
Z µ ¶
h 2 3 (ξh)2 4 (10)
=− Di u(jh) + (ξh)Di u(jh) + Di u(z) ϕ(ξ) dξh
2 2
= −2−1−2n Di2 u(j2−n ) + R
¡ ¢4
with |R| 6 43 h2 kDi4 u(z)k. Here k·k denotes the maximum norm in the hn -neighborhood of
xn,j . Such coefficients, and hence such approximate second derivatives, are directly available
at the hierarchical points (points with odd j). At points with even j (not at the boundary),
we can use
Z
h x
an,2j = − Di2 u(x) ϕ( − 2j) dx
2 h
Z µ ³ ´ ϕ( x − (2j − 1)) + ϕ( x − (2j + 1)) ¶
h 2 x
=− Di u(x) ϕ − 2j − h h
dx (11)
2 2h 2
1
= (an−1,j − an,j−1 − an,j+1 ).
2

The values an−1,j are available (by recursion) from coarser grids, where we have

an−1,j ≈ −21−2n Di2 u(j21−n ) .

This procedure can be used in any of the coordinate directions i. As a result, we find the
nodal representation of the (approximate) second derivative in the i-th direction. If we start
with the H-representation (in every direction) of the function, this procedure delivers the
second derivative (in the i-th direction) in E-representation in the i-th direction and in H-
representation in all other directions. This can simply be converted in either an E- or an
H-representation.
230 P.W. Hemker and F. Sprengel

4. Finite differences on sparse and full grids

As we learn from [13,32], the finite difference operators on sparse grids have a rather complex
structure, including not only finite differences in each direction and linear combinations but
also transformations between the E- and the H-representation. The resulting matrix is ill-
conditioned and even for the Jacobi-preconditioned version rather many BiCGStab iterations
are needed to solve the linear system. On the other hand, finite difference operators on full
grids are well studied. They have a simple structure (7-point stencil for the 3D-Laplacian)
and multigrid routines are at hand. So, the question arises whether and how we can efficiently
use the finite difference operators on the full grids that are contained in the sparse grid in a
multigrid routine for solving the problem on the sparse grid.
In order to answer this question, we analyze the discretized operator (8) as described in
Section 3. For simplicity, we restrict ourselves to the model problem of Poisson’s equation
with homogeneous Dirichlet boundary conditions
−∆u = f in Ω,
(12)
u|δΩ = 0,

on the cube Ω = (0, 1)d and a regular sparse grid. We use the hierarchical representations
of the functions and their second derivatives. For this special case, we give explicit formulas
for the entries of the sparse grid finite difference matrix in hierarchical representation. This
will help us to establish relations between the finite difference discretization on sparse and
full grids. From this, we will finally propose multilevel-type algorithms to solve the linear
system of equations resulting from the finite difference discretization on the sparse grid.

4.1. Characterization of sparse grid and hierarchical points

To introduce the notation for the discussion to follow, it is useful first to investigate the
index sets associated with sparse grid points and hierarchical points, respectively. Consider
a sparse grid Ω+ d
L of level L on Ω = [0, 1] with a finest cell volume h = 2
−L
. Then, every
+
sparse grid point is also a point of a regular grid ΩLe (with step size h in each direction).
Now, we want to characterize the points of the regular grid which belong to the sparse grid,
i.e., we wish to characterize the index set JL with

JL = {` | xLe,` ∈ Ω+
L} .

For this, we use the notation of bit reversing. Let the integer k, satisfying 0 6 k < 2` , have
the binary representation
X`−1
k= ks 2s , ks ∈ {0, 1} . (13)
s=0

Then, the binary representation reflected at position ` − 1 (or in bit reversed order) is given
by1
X`−1
M` (k) = ks 2`−1−s .
s=0

We define the bit reversing of a multi-integer k with respect to a multi-integer ` as M` (k) =

(M`1 (k1 ), . . . , M`d (kd )) and with a scalar m by Mm (k) = Mme (k).
1
For 0 6 k 6 2` , we additionally define M` (2` ) = −1 and M` (−1) = 2` .
 Finite difference operators on dyadic grids 231

For r ∈ Nd0 , we also introduce the set of multi-integers within a dyadic range:

Kr = {k | ki = 2ri −1 , . . . , 2ri − 1, i = 1, . . . , d}

in order to define2 the family of multi-integers with dyadic range up to level L:

[
HL = Kr .
6
|r| L

Then one can easily prove (generalization of [35, 36]) that the mapping of bit reversing
between the two index sets
ML : HL → J L
S
is a bijection. Furthermore, the set |r|=L Kr exactly characterizes the grid points that are
added to the sparse grid of level L − 1 to obtain the sparse grid of level L.
One single set Kr represents (i.e., are bit reversed indices of indices of) the hierarchical
+
points in the subset Ω+r ⊂ ΩL . Thus, the hierarchical points in Ωr are characterized by

{xr,Mr (k) = xLe,ML (k) | k ∈ Kr }.

Remark 1. The same technique can be used starting on a coarsest grid with smaller
mesh width 1/2t0 instead of 1. Then M` has to be defined as a bit reversing w.r.t. ` + t0
instead of `. If we do not start with powers of 2 the situation is more complex but in principle
the same approach can be followed (see [35]).

4.2. Hierarchical representation for finite differences in the univariate case

In equation (11), we gave a recursive expression for the nodal representation of the second
derivative. In this section, we find an explicit expression for the hierarchical coefficients of
the second derivatives found in Section 3.2. With these explicit formulas, we pave the way
for the following sections where explicit formulas for finite differences on full and sparse grids
in hierarchical representation are derived from the univariate expressions by tensor products.
We first develop the hierarchical representation for second order finite differences in the
univariate case. Then, our model problem reads as

−u00 (x) = f for x ∈ (0, 1), u(0) = u(1) = 0.

We interpolate the function u on the equidistant grid with step size hn by a linear spline un .
This spline function has the hierarchical representation
n X
X
un = am,k ϕm,k .
m=1 k odd

Note that, because of the homogeneous boundary conditions, the hierarchical coefficients
vanish on the coarsest level. We approximate the second derivative of a function u by the
second difference
1
∆2,h u(x) = 2 (u(x + h) − 2u(x) + u(x − h)).
h
2
To include multi-integers k with maxi ki = 2` , we denote N` = 2` for ` ∈ N0 and N−1 = −1. and
Kr = {k, ki = Nri −1 , . . . , Nri − 1, i = 1, . . . , d}.
232 P.W. Hemker and F. Sprengel

From this formula for ∆2,hn u, we would obtain the wellknown matrix of second differences
in nodal representation directly. We now want to find a closed form for its hierarchical
representation. In hierarchical points on the finest grid, the hierarchical coefficients represent
the second difference an,k = −2−1−n ∆2,hn u(xn,k ) (cf. Section 3.2). In nonhierarchical points,
we can use the recurrence relation (11). Using bit reversing, we write this now in closed
form. Let xm,k ∈ Ω+ +
n be a hierarchical point (k odd) on the subgrid Ωm with step size hm ,
then Mm (k) ∈ Km and writing out the recursion (11), we obtain
µ n
X ¶
1+n m `−1
−∆2,hn u(xm,k ) = 2 2 am,k − 2 (a`,2`−m k−1 + a`,2`−m k+1 ) . (14)
`=m+1

Now, we interpolate the nodal approximation on level n of the second derivative of the
function u by a linear spline3 ũn and obtain
n −1
2X n
X X
−ũn = −∆2,hn u(xn,k ) ϕn,k = bm,k ϕm,k . (15)
k=1 m=1 k odd

This defines the coefficients bm,k which clearly satisfy

1
bm,k = (∆2,hn u(xm,k−1 ) − 2∆2,hn u(xm,k ) + ∆2,hn u(xm,k+1 ))
2
for hierarchical points xm,k . The nonhierarchical points xm,k−1 and xm,k+1 appear in this
formula, too. If we want to apply formula (14), we need to characterize these points (non-
hierarchical on level m) as hierarchical points on some coarser level. For this, we need some
further notation.
Let k have the binary representation (13). Now we need to characterize the indices m−
and m+ for which Mm (k − 1) ∈ Km− and Mm (k + 1) ∈ Km+ , where m± denotes the level on
which xm,k±1 is hierarchical. We see that
(
0 for k = 1,
m− =
max {m − 1 − s, ks = 1} otherwise,
s=1,...,m−1

and (
0 for k = 2m − 1,
m+ =
max {m − 1 − s, ks = 0} otherwise.
s=1,...,m−1

Thus, if the points xm,k are not next to the boundary, i.e., for4 m± 6= 0,
µ
1+n
−∆2,hn u(xm,k±1 ) = 2 2m± am± ,2m± −m (k±1)
n
X ¶
`−1
− 2 (a`,2`−m (k±1)−1 + a`,2`−m (k±1)+1 ) .
`=m± +1

3
For notational convenience, we set in the end points ∆2,hn u(0) = ∆2,hn u(1) = 0.
4
For points next to the boundary, i.e., for m± = 0, we use the notation from the last three footnotes and
set ∆2,hn u(xm,k±1 ) = 0.
 Finite difference operators on dyadic grids 233

Collecting all terms, we get the hierarchical coefficient

µ Xn
1+n m
bm,k = 2 2 am,k − 2`−1 (a`,2`−m k−1 + a`,2`−m k+1 )
`=m+1
n
X
m− −1
−2 am− ,2m− −m (k−1) + 2`−2 (a`,2`−m (k−1)−1 + a`,2`−m (k−1)+1 ) (16)
`=m− +1
Xn ¶
m+ −1 `−2
−2 am+ ,2m+ −m (k+1) + 2 (a`,2`−m (k+1)−1 + a`,2`−m (k+1)+1 )
`=m+ +1

with modifications for m± = 0. Note that this expression depends on m, k, and n, i.e., it
depends not only on the point xm,k but also on the finest grid chosen.
We can write the coefficients am,k and bm,k in vector notation and describe the trans-
formation process by a matrix. For this, we combine hierarchical coefficients in a vector
as ¡ ¢
v n = vm,k 16m6n,Mm (k)∈Km (17)
in order to define the vectors an and bn , and nodal coefficients simply as
¡ ¢T
un = u(xn,1 ), . . . , u(xn,2n −1 ) .

In this way, we define the matrices An and H n by

An an = bn and H n un = an .

By construction, H −1n An H n is the usual matrix of second differences in nodal representation

and hence symmetric positive definite. Thus, the matrix An of finite differences in terms of
hierarchical coefficients has the same, i.e., only real positive eigenvalues. However, it is not
symmetric and cannot be symmetric positive definite.
Solving the discretized system in its hierarchical form now comes down to solving
¡ ¢T
An an = H n f (xn,1 ), . . . , f (xn,2n −1 ) .

4.3. Hierarchical representation for finite differences on the sparse grid

To gain some more insight into the structure of the sparse grid FD operator, we want to
obtain an explicit expression for the hierarchical coefficients of the discrete Laplacian. That
is, we look for an expression for the elements of the FD matrix AL . To this end, we investigate
the discrete second derivative in the x1 -direction first.
Let L be the (highest) level of the sparse grid. Let uL be given in H-representation:
X X
uL = am,` ϕm,` .
6
|m| L |||`||| odd

Because of the boundary conditions, here and in the sequel the notation |m| 6 L means:
m > 0, |m| 6 L. (The coefficients responsible for the boundary are zero.) Let ũL be the
approximation of the resulting function ∆u on the sparse grid Ω+
L with
X X
−ũL = b∆
m,` ϕm,` .
6
|m| L |||`||| odd
234 P.W. Hemker and F. Sprengel

We denote the approximation of the second derivative in the xν -direction on the sparse grid
Ω+
L by X X (ν)
(ν)
−ũL = bm,` ϕm,` .
6
|m| L |||`||| odd

Then, the hierarchical coefficients obviously fulfill

(1) (2) (d)
b∆
m,` = bm,` + bm,` + · · · + bm,` . (18)

Let xm,k ∈ Ω+ + +
L be a hierarchical point on the grid Ωm . The full grid Ωj ⊂ ΩL , which
is the finest in the x1 -direction such that xm,k ∈ Ω+ j , is characterized by the multi-index
j = (m1 + L − |m|, m2 , . . . , md ). Then the hierarchical coefficient of the second difference in
the x1 -direction at point xm,k is

(1) 1¡ ¢
bm,k = −αm,k + αm,k−e1 + αm,k+e1 , (19)
2
where αm,` denotes the coefficient of the approximation of the second derivative in the x1 -
direction (E-representation in the x1 -direction, H-representation in all other directions) at
the point xm,` .
By construction Mm (k) ∈ Km . Now using the results of the previous section, we obtain
from the finite difference operator in the x1 -direction, see (14),
µ m1 +L−|m|
X
1+2m1 +L−|m|
¡
αm,k = −2 am,k − 2`−m1 −1 a(`,m2 ,...,md ),(2`−m1 k1 −1,k2 ,...,kd )
`=m1 +1
¶
¢
+ a(`,m2 ,...,md ),(2`−m1 k1 +1,k2 ,...,kd ) .

In the same way as in the previous section, we denote by m1− and m1+ the indices for which
Mm1 (k1 − 1) ∈ Km1− and Mm1 (k1 + 1) ∈ Km1+ (i.e., m1− and m1− are the x1 -levels on which
xm,k are hierarchical points in the x1 -direction). Now we can characterize the remaining
terms in (19) and obtain
µ m1 +L−|m|
X
(1) 1+2m1 +L−|m|
¡
bm,k =2 · am,k − 2`−m1 −1 a(`,m2 ,...,md ),(2`−m1 k1 −1,k2 ,...,kd )
`=m1 +1
¢
+ a(`,m2 ,...,md ),(2`−m1 k1 +1,k2 ,...,kd )
− 2m1+ −m1 −1 a(m1+ ,m2 ,...,md ),(2m1+ −m1 (k1 +1),k2 ,...,kd )
m1 +L−|m|
X ¡
+ 2`−m1 −2 a(`,m2 ,...,md ),(2`−m1 (k1 +1)−1,k2 ,...,kd )
`=m1+ +1
¢ (20)
+ a(`,m2 ,...,md ),(2`−m1 (k1 +1)+1,k2 ,...,kd )
− 2m1− −m1 −1 a(m1− ,m2 ,...,md ),(2m1− −m1 (k1 −1),k2 ,...,kd )
m1 +L−|m|
X ¡
+ 2`−m1 −2 a(`,m2 ,...,md ),(2`−m1 (k1 −1)−1,k2 ,...,kd )
`=m1− +1
¶
¢
+ a(`,m2 ,...,md ),(2`−m1 (k1 −1)+1,k2 ,...,kd )
 Finite difference operators on dyadic grids 235

(with obvious modifications for points next to the boundary, where m1− = 0 or m1+ = 0).
Here, the finest grid in the x1 -direction depends on the other coordinates of the evaluation
(1)
point. Therefore, bm,k depends on m, k and the highest level L.
(ν) ∆ (ν)
Now, we write the hierarchical coefficients bm,k , b∆
m,k , and am,k in vector form as bL , bL ,
and aL , respectively, using the d-dimensional version of (17):
¡ ¢
v L = vm,k |m|6L,Mm (k)∈Km .

(ν)
We define the matrices AL and AL by
(ν) (ν)
AL aL = b∆
L and AL aL = bL , (21)
(1) (2) (d)
respectively. Then, obviously AL = AL + AL + · · · + AL .

Remark 2. Note that the formulas (16), (20) and the corresponding formula on full
grids can also be used in the case of adaptive grids and deliver another representation of the
same finite differences as in [32]. The sums run over the hierarchical points from the coarser
to the finer grids. Leaving out points on finer grids (if the hierarchical coefficients are set to
zero) only shortens the sums. The restriction to regular sparse grids is only for simplicity of
presentation.

4.4. Hierarchical representation for finite differences on full grids

The finite difference discretization of (12) on a full grid is wellknown (the usual 7-point stencil
in 3D). Here we rewrite its matrix in the hierarchical representation in order to unfold the
relation between the finite difference operators on full and sparse grids, with the aim to
find proper multigrid algorithms. Let a discrete function on Ω+ n , say un , be given in its
hierarchical representation
X X
un = am,k ϕm,k .
6
0<m n |||k||| odd

Let ũ∆
n denote the approximation of ∆u, similar to (15) in the one-dimensional case in
Section 4.2, X X
−ũ∆n = b̃∆
m,k ϕm,k .
6
0<m n |||k||| odd

(ν)
Assume that b̃m,k are hierarchical coefficients in the approximation of the second derivative
(ν) +
in the xν -direction ũn . The coefficients corresponding to the subgrid Ω+
n ⊂ ΩL are written
in vector-form as ¡ ¢
v n = vm,k m6n,Mm (k)∈Km .
∆ (1) (2) (d)
Similar to (18) in the previous section, we have b̃n = b̃n + b̃n + · · · + b̃n , and we write
(ν) ∆
A(ν)
n an = b̃n and An an = b̃n .

Obviously,
An = A(1) (2) (d)
n + An + . . . + An . (22)
236 P.W. Hemker and F. Sprengel

S © ª
We introduce a one-dimensional index set by In = nm=1 (m, k) | Mm (k) ∈ Km and the
corresponding identity matrix by
¡ ¢#(In )
I n = δµ,ν µ,ν=1 .
© ª
Then, for the d-dimensional case (m, k) | Mm (k) ∈ Km = In1 × In2 × · · · × Ind and we
find the matrix of second differences in the ν-th direction as the Kronecker product

A(ν)
n = I n1 ⊗ · · · ⊗ I nν−1 ⊗ Anν ⊗ I nν+1 ⊗ · · · ⊗ I nd .

where Anν are the same as in Section 4.2. As a Kronecker sum of the matrices Anν , the
matrix An also has only positive eigenvalues. This is not surprising because An has the
same eigenvalues as the finite difference matrix for the nodal representation on the full grid
which is symmetric positive definite.

4.5. Relations between finite differences on full and sparse grids

In this section, we establish relations between the discrete Laplacian on full and on sparse
grids. If one considers the Galerkin or FEM approach for discretization on sparse grids,
using the hierarchical basis or the generating system of nodal bases functions as test and
trial functions, one obtains simple relations between the stiffness matrices for the sparse
+
and the full grids. The basis functions for a full grid Ω+ n ⊂ ΩL simply form a subset of the
hierarchical basis/generating system for the sparse grid. So, by construction it can be written
as a Galerkin product using the matrix for the sparse grid. As a result, one can immediately
write down multiplicative subspace correction algorithms like in [15] and interpret them as
block iteration methods.
In the finite difference case things are different. In contrast to the sparse grid, on a full
grid, for a constant coefficient differential equation the evaluation of the finite difference at
a certain point does not depend on the location of the point. For this reason, the finite
difference matrix of a full grid Ω+ n cannot be written as a Galerkin product including the
finite difference matrix of the sparse grid Ω+ L of level L > |n|. This will become immediately
clear from equation (24), below.
In what follows, we discuss the kind of relation between the finite difference operators for
the Laplacian discretized on a full and on a sparse grid, and we propose a solution algorithm
resulting from this.

4.5.1. Full grids as subgrids of the sparse grid. Assume |n| 6 L. We define a
prolongation P L,n : Vn → VL by wL = P L,n v n with
(
vm,k for m 6 n,
wm,k =
0 otherwise,

where vm,k and wm,k are hierarchical coefficients. This represents linear interpolation on
(ν)
the sparse grid. Further, we need a direction-dependent restriction Rn,L : VL → Vn by
(ν)
wn = Rn,L v L with
wm,k = 2nν −mν −L+|m| vm,k for m 6 n . (23)
The scaling factor in the above definition arises from the difference in scaling between (16)
and (20). Comparing these two formulas and keeping (21) and (22) in mind, we see that,
 Finite difference operators on dyadic grids 237

(ν)
with these P L,n and Rn,L , we can write the matrix of finite differences for the Laplacian
on the full grid as a sum of Galerkin products using the direction (and point) dependent
weighted finite differences for the second derivatives in each direction:
¡ (1) (1) (2) (2) (d) (d) ¢
An = Rn,L AL + Rn,L AL + · · · + Rn,L AL P L,n . (24)

Unfortunately, this representation is not in a Galerkin form suited for multigrid-type algo-
rithms. We can see from (24) that the finite difference operators on full grids themselves
cannot serve as good approximations for the finite difference operator on the sparse grid.

4.5.2. Full grids as subgrids of the sparse grid — An alternative approach.

(ν)
In the previous paragraph, we scaled the different parts, AL , of the discrete Laplacian for
(ν)
the sparse grid for the different directions. On the other hand, we can scale An , i.e., the
directional parts of the matrix responsible for the full grid. With the scaling matrices
¡ mν −nν ¢
M (ν)
n = diag 2 m6n,Mm (k)∈Km
¡ ¡ ¢ ¢
= I n1 ⊗ · · · ⊗ I nν−1 ⊗ diag 2mν −nν (mν ,kν )∈In ⊗ I nν+1 ⊗ · · · ⊗ I nd
ν

and (22), we define

Ãn = M (1) (1) (2) (2) (d) (d)

n An + M n An + · · · + M n An . (25)

This matrix is again a Kronecker sum of matrices with positive eigenvalues and so it has
only positive eigenvalues. Introducing another, differently scaled, restriction that is now
independent of the direction ν, R̃n,L : VL → Vn , defined by wn = R̃n,L v L with

wm,k = 2−L+|m| vm,k for m 6 n ,

the new matrix can be written as the Galerkin product

Ãn = R̃n,L AL P L,n .

This Galerkin relation may lead to the following solution algorithm which uses the ma-
trices Ãn . Denote the vector of hierarchical coefficients of the right-hand side by f L , and let
the approximation aL of the hierarchical coefficients of uL be given. Then the next iteration
step is
for all |n| = k
−1
do aL := aL + P L,n Ãn R̃n,L (f L − AL aL ) (26)
enddo
This defect correction algorithm of multiplicative type corresponds to a block iterative solver
for the discrete Laplacian on the sparse grid with overlapping blocks. Seen as such a block
iterative solver this is to some extent similar to the FEM case (see [15]), but in the case
of finite differences the blocks (containing point dependent scaling) cannot really be seen
as a finite difference discretization on a certain full grid. Nevertheless, in the special case
of the Laplacian it is still possible to give a relation (25) between the blocks and the finite
difference discretizations.
To improve the convergence rate of the above algorithm we might think of using more
than one level and/or more than one iteration per level which leads to
238 P.W. Hemker and F. Sprengel

cycles with L = 6,...,9, l = L cycles with L =6,...,9, l = 3

1e1 1e1
L=6 L=6
5e0 L=7 5e0 L=7
L=8 L=8
L=9 L=9

1e0 1e0

5e–1 5e–1

Residual Residual

1e–1 1e–1

5e–2 5e–2

1e–2 1e–2

5e–3 5e–3

2 4 6 8 10 12 14 16 18 2 4 6 8 10 12 14 16 18
Number of cycles Number of cycles

Figure 1. Left: Convergence of Algorithm (26) for the levels L = 6, . . . , 9. Right: Convergence of Algorithm
(27) for the levels L = 6, . . . , 9, ` = 3 with ν = 1.

for k = ` to L
do for i = 1 to ν
do for all |n| = k
−1
do aL := aL + P L,n Ãn R̃n,L (f L − AL aL ) (27)
enddo
enddo
enddo

with a fixed lower level ` > d and a number ν of iterations per level.

Example: We apply the Algorithms (26) and (27) to the following 3D-problem. Solve
(12) with the right-hand side

³Y
3 3
Y ´
2
f (x) = −3π sin πxi + 64 sin 8πxi
i=1 i=1

(0)
and starting from the zero function uL ≡ 0. We obtain the convergence behavior shown
in Fig. 1. We see that we get better convergence if we include also lower levels (right). In
both cases, the speed of convergence slows down with growing level. Approximately, the
reduction factor gets worse with L2 , the square of the highest level.
If we apply Algorithms (26) and (27) to a similar 2D-examples, we obtain reduction
factors that grow linearly with L. This is comparable to [32] where the iteration number for
the Jacobi-preconditioned BiCGStab for similar 2D-problems grow at least linear in L. In a
recent paper [14], a further reduction of the necessary BiCGStab iterations could be achieved
by wavelet preconditioning, but also in this case the growth of the iteration numbers seems
to be linear in L.

Discussion
Since the efficient O(N ) evaluation of finite element stiffness matrices for variable coefficient
equations on sparse grids (d > 2) is impossible, and FE residual computation requires
unattractive modifications of the FE discretization [9] (cf. [40]), one might be tempted to use
 Finite difference operators on dyadic grids 239

finite differences (FD) instead. With the aim to construct an efficient multigrid procedure,
we analyze the structure of the FD discrete Laplacian, and Galerkin relations are established.
In particular, we study the relation between the FD operators on full and sparse grids. The
rather complex relations (24) and (25) for this simple special case show that it is unlikely
that – for the general case – a natural and efficient multilevel preconditioner can be found
that shows spectral equivalence.
Thus, so far all attempts to construct O(N ) solution methods for finite differences on
ASG failed, and the structure of the discretization as analyzed in this paper gives little hope.
It is our conclusion that, for the solution of the general equation (8) in three dimensions, it
is probably better to abandon the generality of sparse grids and to adopt condition L. Then,
representations on coarse grids are simply obtained by restriction from a unique finest grid,
the application of FE gives no particular problems and the usual Galerkin relations hold.
Optimal efficiency can be obtained by exploitation of all C k -regularity and by the appli-
cation of hp-adaptive methods [2] (with dyadic refinement) and multigrid semi-coarsening
techniques.

References
[1] R. Balder, Adaptive Verfahren für elliptische und parabolische Differentialgleichungen auf dünnen Git-
tern, Ph.D. thesis, Technische Universität München, 1994.

[2] T. J. Barth and H. Deconinck, High-Order Methods for Computational Physics, vol. 9 of Lecture Notes
in Computational Science and Engineering, Springer Verlag, Harlow, 1999.

[3] V. Barthelmann, E. Novak, and K. Ritter, High dimensional polynomial interpolation on sparse grids,
Adv. Comput. Math., 12 (1999), pp. 273–288.

[4] H.-J. Bungartz and T. Dornseifer, Sparse grids: Recent developments for elliptic partial differential equa-
tions, in: Multigrid Methods V (W. Hackbusch and G. Wittum, eds.), Lecture Notes in Computational
Science and Engineering 3, Springer, Berlin, 1998.

[5] H.-J. Bungartz and M. Griebel, A note on the complexity of solving Poisson’s equation for spaces of
bounded mixed smoothness, J. Complexity, 15 (1999), pp. 167–199.

[6] H.-J. Bungartz, M. Griebel, D. Röschke, and C. Zenger, Pointwise convergence of the combination
technique for the Laplace equation, East-West J. Numer. Math., 2 (1994), pp. 21–45.

[7] F.-J. Delvos and W. Schempp, Boolean Methods in Interpolation and Approximation, Pitman Research
Notes in Mathematics Series 230, Longman Scientific & Technical, Harlow, 1989.

[8] R. A. DeVore, S. V. Konyagin, and V. N. Temlyakov, Hyperbolic wavelet approximation, Constr. Ap-
prox., 14 (1998), pp. 1–26.

[9] T. Dornseifer, Diskretisierung allgemeiner elliptischer Differentialgleichungen in krummlinigen Koordi-

natensystemen auf dünnen Gittern, Ph.D. thesis, Technische Universität München, 1997.

[10] T. Gerstner and M. Griebel, Numerical integration using sparse grids, Numer. Algorithms, 18 (1998),
pp. 209–232.

[11] M. Griebel, A parallelizable and vectorizable multilevel algorithm on sparse grids, in: Parallel algorithms
for partial differential equations, Notes Numer. Fluid Mech. 31, Vieweg, Wiesbaden, 1991, pp. 94–100.

[12] M. Griebel, Multilevelmethoden als Iterationsverfahren über Erzeugendensystemen, Teubner Skripten

zur Numerik, Teubner, Stuttgart, 1994.

[13] M. Griebel, Adaptive sparse grid multilevel methods for elliptic PDEs based on finite differences, Com-
puting, 61 (1998), pp. 151–180.
240 P.W. Hemker and F. Sprengel

[14] M. Griebel and F. Koster, Adaptive wavelet solvers for the unsteady incompressible Navier–Stokes equa-
tions, in: Advances in Mathematical Fluid Mechanics, (J. Malek, J. Necas, and M. Rokyta, eds.),
Lecture Notes of the Sixth International School Mathematical Theory in Fluid Mechanics, Paseky,
Czech Republic, September 1999. Springer Verlag, 2000.
[15] M. Griebel and P. Oswald, On the abstract theory of additive and multiplicative Schwarz algorithms,
Numer. Math., 70 (1995), pp. 163–180.
[16] M. Griebel and T. Schiekofer, An adaptive sparse grid Navier–Stokes solver in 3D based on finite
differences, in: Proceedings ENUMATH97 (H.G. Bock, G. Kanschat, R. Rannacher, F. Brezzi, R.
Glowinski, Y.A. Kuznetsov, and J. Periaux, eds.), Heidelberg, 1999. World Scientific Publishing.
[17] M. Griebel, M. Schneider, and C. Zenger, A combination technique for the solution of sparse grid
problems, in: Iterative Methods in Linear Algebra (R. Beauwens and P. de Groen, eds.), Elsevier,
IMACS, Amsterdam, 1992, pp. 263–281.
[18] P. W. Hemker, Sparse–grid finite–volume multigrid for 3D–problems, Adv. Comput. Math., 4 (1995),
pp. 83–110.
[19] P. W. Hemker and P. de Zeeuw, BASIS3, a data structure for 3-dimensional sparse grids, in: Eu-
ler and Navier-Stokes Solvers Using Multidimensional Upwinds Schemes and Multigrid Acceleration
(H. Deconinck, ed.), Notes Numer. Fluid Mech. 57, Vieweg, Wiesbaden, 1997, pp. 443–484.
[20] P. W. Hemker and C. Pflaum, Approximation on partially ordered sets of regular grids, Appl. Numer.
Math., 25 (1997), pp. 55–87.
[21] P. W. Hemker and F. Sprengel, Experience with the solution of a finite difference discretization on
sparse grids, in: Numerical Analysis and Its Applications (L. Vulkov, J. Wasniewski, and P. Yalamov,
eds.), Springer Verlag, Berlin, 2001, pp. 402–413, Proc. of the Conference “Numerical Analysis and its
Applications”, 2000, Rousse.
[22] N. M. Korobov, Approximate calculation of multiple integrals with the aid of methods in the theory of
numbers, Dokl. Akad. Nauk SSSR, 115 (1957), pp. 1062–1065, in Russian.
[23] F. Kupka, Sparse Grid Spectral Methods for the Numerical Solution of Partial Differential Equations
with Periodic Boundary Conditions, Ph.D. thesis, Universität Wien, 1997.
[24] W. Mulder, A high-resolution Euler solver based on multigrid semi-coarsening, and defect correction,
Journal of Computational Physics, 100 (1992), pp. 91–104.
[25] W. A. Mulder, A new multigrid approach to convection problems, Journal of Computational Physics,
(1989), pp. 303–323.
[26] J. Noordmans and P. W. Hemker, Convergence results for 3D sparse grid approaches, Numer. Linear
Algebra Appl., 5 (1999), pp. 363–376.
[27] J. Noordmans and P. W. Hemker, Application of an adaptive sparse-grid technique to a model singular
perturbation problem, Computing, 65 (2000) pp. 357–378.
[28] E. Novak and K. Ritter, High dimensional integration of smooth functions over cubes, Numer. Math.,
75 (1996), pp. 79–98.
[29] P. Oswald, Multilevel Finite Element Approximation. Theory and Applications, Teubner Skripten zur
Numerik, Teubner, Stuttgart, 1994.
[30] C. Pflaum, Diskretisierung elliptischer Differentialgleichungen mit dünnen Gittern, Ph.D. thesis, Tech-
nische Universität München, 1996.
[31] C. Pflaum, A multilevel algorithm for the solution of second order elliptic differential equations on sparse
grids, Numer. Math., 79 (1998), pp. 141–155.
[32] T. Schiekofer, Die Methode der Finiten Differenzen auf dünnen Gittern zur Lösung elliptischer und
parabolischer partieller Differentialgleichungen, Ph.D. thesis, Universität Bonn, 1998.
[33] T. Schiekofer and G. Zumbusch, Software concepts of a sparse grid finite difference code, in: Concepts
of Numerical Software (W. Hackbusch and G. Wittum, eds.), Notes in Numerical Fluid Mechanics,
Vieweg, Braunschweig, 1998, .
 Finite difference operators on dyadic grids 241

[34] S. A. Smolyak, Quadrature and interpolation formulas of the classes Wsa and Esa , Dokl. Akad. Nauk
SSSR, 131 (1960), pp. 1028–1031, in Russian, Engl. Transl.: Soviet Math. Dokl. 1:384–387, 1963.
[35] F. Sprengel, Interpolation and Wavelet Decomposition of Multivariate Periodic Functions, Ph.D. thesis,
University of Rostock, 1997.
[36] F. Sprengel, Interpolation and wavelets on sparse Gauß–Chebyshev grids, in: Multivariate Approxi-
mation, Recent Trends and Results (W. Haußmann, K. Jetter, and M. Reimer, eds.), Mathematical
Research, Vo. 101, Akademie–Verlag, Berlin, 1997, pp. 269–286.
[37] F. Sprengel, Interpolation of functions from Besov–type spaces on Gauß–Chebyshev grids, J. Complexity,
16 (2000), pp. 507-523.
[38] C. Teitzel, M. Hopf, R. Grosso, and T. Ertl, Volume visualization on sparse grids, Computing and
Visualization in Science, 2 (1999), pp. 47–59.
[39] V. N. Temlyakov, Approximation of Periodic Functions, Nova Science, New York, 1993.
[40] E. E. P. Veldman and R. W. C. P. Verstappen, Higher–order discretization methods for CFD, Nieuw
Archief voor Wiskunde, Ser. IV Vol.17 (1999), pp. 195–204.
[41] G. W. Wasilkowski and H. Woźniakowski, Explicit cost bounds of algorithms for multivariate tensor
product problems, J. Complexity, 11 (1995), pp. 1–56.
[42] C. Zenger, Sparse grids, in: Parallel Algorithms for Partial Differential Equations (W. Hackbusch, ed.),
Notes on Numerical Fluid Mechanics 31, Vieweg, Braunschweig, 1991, pp. 297–301.

Received 20 Feb. 2001

Revised 20 Jul. 2001