# The numerical performance of wavelets for PDEs: the multi-scale finite element

Mark Christon, David Roach

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mark Christon, David Roach. The numerical performance of wavelets for PDEs: the multi-scale finite element. Computational Mechanics, 2000, 25, pp.230-244. 10.1007/s004660050472. hal-01333698

HAL Id: hal-01333698

https://hal.science/hal-01333698

Submitted on 18 Jun 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# The numerical performance of wavelets for PDEs: the multi-scale finite element

**M. A. Christon**

Computational Physics R&D Department,
Sandia National Laboratories,
M/S 0819, P.O. Box 5800 Albuquerque,
New Mexico 87185-0819, USA

**D. W. Roach**

Mathematics Department, University of Georgia, USA

Sandia is a multiprogram laboratory operated by Sandia Corporation for the United States Department of Energy under Contract DE-AC04-94AL85000.

**Abstract** The research summarized in this paper is part of a multi-year effort focused on evaluating the viability of wavelet bases for the solution of partial differential equations. The primary objective for this work has been to establish a foundation for hierarchical/wavelet simulation methods based upon numerical performance, computational efficiency, and the ability to exploit the hierarchical adaptive nature of wavelets. This work has demonstrated that hierarchical bases can be effective for problems with a dominant elliptic character. However, the strict enforcement of orthogonality in the usual $L^2$ sense is less desirable than orthogonality in the energy norm. This conclusion has led to the development of a multi-scale linear finite element based on a hierarchical change-of-basis. This work considers the numerical and computational performance of the hierarchical Schauder basis in a Galerkin context. A unique row-column lumping procedure is developed with multi-scale solution strategies for 1-D and 2-D elliptic partial differential equations.

## 1

### Introduction

Wavelets are a relatively new mathematical tool that distinct data, functions, and differential operators into components with an associated resolution that is matched to the scale of each component. The name “wavelet” or “ondelette” was coined in the early 1980s by French researchers Morlet, Arens, Fourgeua, Giard and Grossman [1, 3, 4, 9]. However, functions with the attributes of wavelets have been known for almost 100 years. Meyer [7] points out that there are seven primary origins for wavelets that date from around 1930 with the Haar wavelet dating back to 1909. However, the literature from this era does not explicitly use the term “wavelet”.

In recent years, the application of wavelet bases to the solution of partial differential equations (PDEs) has evolved to the point where there are a number of com-

peting formulations that include, but are not limited to wavelet-Galerkin, wavelet-collocation, and reproducing kernel methods. Despite the growing number of wavelet-based formulations and solution algorithms for PDEs, the field is still very new, and many technical issues remain. That is to say, there is promise in the approach, but there is a clear need for fundamental research that characterizes the numerical and computational performance of wavelets for the solution of partial differential equations.

A complete historical perspective of wavelet based approaches to solving PDEs is beyond the scope of this paper. However, a survey of the literature on wavelet based methods for solving PDEs may be found in Christon et al. [2]. Over the past 8 years, there has been a great deal of attention paid to the wavelet-Galerkin method because of the generality that the Galerkin method provides, and the ease with which alternative bases may be implemented and tested. In addition, the trend in recent years has been away from traditional Daubachies wavelets with $L^2$ orthogonality and towards wavelets that are constructed with weaker orthogonality properties yielding biorthogonal and semi-orthogonal wavelets (see for example Dahmen and Micchelli [4]). This approach has led to the idea that wavelets may be most useful for solving PDEs if they are used to simply “complete” a subspace with the concomitant relationship between grid scales. A related approach that combines a spectral method with biorthogonal wavelet bases adapted to non-square domains and higher spatial dimensions may be found in Canuto et al. [1].

In this work, orthogonality with respect to the energy norm is used with a nodal change-of-basis to develop multi-scale finite elements that make use of a hierarchical Schauder basis. In 1-D, the uniform stability of the Schauder basis in $H^1$ yields uniformly bounded condition numbers, i.e., independent of refinement level (grid resolution), for elliptic operators in a Galerkin form. This property was previously described by Zienkiewicz [12] as the “diagonality” feature for one dimensional hierarchical shape functions. However, the relationship between uniform stability, a nodal change-of-basis and the wavelet transform was not considered. In the ensuing discussion, the theoretical issues surrounding the hierarchical Schauder basis are presented, and the details of a uniformly stable basis are made concrete. The use of the wavelet transform in conjunction with the assembly of the mass and stiffness operators for a hierarchical basis is presented in the framework of 1-D and 2-D multi-scale finite elements. In addition, a comparison between the hierarchical Schauder basis and the linear finite element
basis is presented. The application of an ad-hoc row-column lumping procedure for the hierarchical basis is outlined, and its effectiveness for solving elliptic boundary value problems in one and two-dimensions is demonstrated.

## 2 The hierarchical Schauder basis

In this section, the formulation issues associated with the hierarchical Schauder basis are outlined. For this development, the boundary value problem under consideration is

$$ -\varepsilon \nabla^2 u + \gamma u = f \quad \text{on } \Omega, \qquad (1) $$

with homogeneous Dirichlet boundary conditions

$$ u = 0 \quad \text{on } \Gamma. \qquad (2) $$

The weak form of the problem, after introducing the finite dimensional subspace, $\mathcal{V}_k$, is

$$ \int_{\Omega} \varepsilon \nabla u_k \cdot \nabla \mathbf{\Phi}_k \, d\Gamma + \int_{\Omega} \gamma u_k \mathbf{\Phi}_k \, d\Gamma = \int_{\Gamma} f \mathbf{\Phi}_k \, d\Gamma. \quad (3) $$

Here, $u_k \in \mathcal{H}_0^1$, and $\mathbf{\Phi}_k$ is a finite dimensional basis for the subspace $\mathcal{V}_k$ at scale $k$. In matrix form, this may be written as

$$ [\mathbf{M}_k + \mathbf{K}_k]\{u_k\} = \{F_k\} \qquad (4) $$

where $\mathbf{M}_k$ is the mass matrix, and $\mathbf{K}_k$ is the stiffness matrix associated with scale $k$. Alternatively, in a bilinear form, $a(\mathbf{\Phi}_k, \mathbf{\Phi}_k)u_k = F(\mathbf{\Phi}_k)$.

Now, consider an alternative basis for $\mathcal{V}_k$,

$$ \mathbf{\Psi}_k = \mathbf{W}_k^\mathrm{T} \mathbf{\Phi}_k, \qquad (5) $$

where $\mathbf{W}_k$ is a nonsingular $N \times N$ matrix. From this transformation, $u_k = \Delta u^\mathrm{T} \mathbf{\Psi}$ may also be found by solving

$$ a(\mathbf{\Psi}_k, \mathbf{\Psi}_k) \Delta u = F(\mathbf{\Psi}_k) \qquad (6) $$

for the multi-scale representation of the field, $\Delta u$. The linear system, $A_k^\Psi$, that derives from the change of basis may also be viewed as the linear system obtained by preconditioning $A_k^\Phi = [\mathbf{M}_k + \mathbf{K}_k]$ with $\mathbf{W}_k$, i.e.,

$$ A_k^\Psi = \mathbf{W}_k^\mathrm{T} \mathbf{A}_k^\Phi \mathbf{W}_k. \qquad (7) $$

Note that the superscript $\Psi$ ($\Phi$) indicates that $\mathbf{\Psi}_k$ ($\mathbf{\Phi}_k$) is used as both the test and trial function in obtaining the weak form.

The matrix $\mathbf{W}_k$ is a wavelet transform and provides the mechanism for decomposing a field into multiple scales. In order to understand the role of the wavelet transform, consider that

$$ \mathcal{V}_0 \subset \mathcal{V}_1 \subset \cdots \subset \mathcal{V}_k \subset \cdots \qquad (8) $$

is a one-sided sequence of nested finite-dimensional subspaces of $\mathcal{H}$, a Hilbert space, such that $\bigcup_{k=0}^\infty \mathcal{V}_k = \mathcal{H}$. Now, define the “coarse-grid” wavelet subspace, $\mathcal{W}_0 = \mathcal{V}_0$ and, for $k \ge 1$, choose $\mathcal{W}_k$ so that completion of subspace $\mathcal{V}_{k-1}$ may be achieved via a direct sum,

$$ \mathcal{V}_k = \mathcal{V}_{k-1} \oplus \mathcal{W}_k. \qquad (9) $$

Here $\oplus$ denotes a direct sum, but not an orthogonal direct sum.

Beginning with Eq. (5), the wavelet transform may be written recursively in terms of a two-scale transform as

$$ \mathbf{W}_k = T_k \begin{bmatrix} T_{k-1} & 0 \\ 0 & I_{k-1} \end{bmatrix} \cdots \begin{bmatrix} T_1 & 0 \\ 0 & I_1 \end{bmatrix} \qquad (10) $$

where $I_k$ is the $N_\Psi \times N_\Psi$ identity matrix with $N_\Psi = \dim(\mathcal{W}_k)$. Here, $T_k$ is a two-scale transform such that

$$ \left\{ \begin{matrix} \mathbf{\Phi}_{k-1} \\ \mathbf{\Psi}_k \end{matrix} \right\} = T_k^\mathrm{T} \mathbf{\Phi}_k. \qquad (11) $$

As an aside, one very important aspect of this formulation is that it does not require the decomposition matrix $\mathbf{W}_k^{-1}$. This is significant because it permits the relaxation of strict orthogonality in the selection of the wavelet bases.

One choice for $\mathcal{W}_k$ that satisfies Eq. (9) is the hierarchical Schauder basis [3, 10–13]. This basis consists of the elements of the linear finite element basis at multiple scales where the relationship between $\mathbf{\Psi}_k$ and $\mathbf{\Phi}_k$ is simply

$$ \psi_k^j = \phi_k^{2j-1}, \quad j = 1, \dots, \mathrm{Nnp}_k, \qquad (12) $$

where $\mathbf{\Phi}_k$ is the linear finite element basis. Here the superscript $j$ indicates the node number, $k$ indicates the scale, and $\mathrm{Nnp}_k$ is the number of node points at scale $k$. The Schauder basis is illustrated in Fig. 1.

With the relationship between basis elements in the Schauder basis given by Eq. (12), the construction of the two-scale transformation may be written in terms of the grid spacing. In order to define the wavelet transform, let $\mathbf{H}_k$ be the $(\mathrm{Nnp}_k - 1) \times (\mathrm{Nnp}_{k-1} - 1)$ matrix $\mathbf{H}_k = (h_k^{j,j})_{j,j}$ and let $\mathbf{G}_k$ be the $(\mathrm{Nnp}_k - 1) \times (\mathrm{Nnp}_{k-1})$ matrix

<figure><img src="image_3.png" /><figcaption>Fig. 1. The hierarchical Schauder basis for scale $k=2$ with 8 elements and 9 nodes showing a the equivalent finite element basis, b the composite basis, c the basis functions for Scale-2, d the basis functions for Scale-1, and e the basis functions for Scale-0</figcaption></figure>
$G_k = (g_k^{j'j})_{j'j}$. Then, the two-scale transformation for the Schauder basis is given by

$$T_k = [\mathbf{H}_k | \mathbf{G}_k] . \qquad (13)$$

The function values $h_k^{j'j}$ and $g_k^{j'j}$ for $\phi_{k-1}^j$ and $\psi_k^j$ respectively at the node points $x_k^{j'}$ are given by

$$h_k^{j'j} = \begin{cases} \frac{\Lambda_k^{j'}}{\Delta_{k-1}^j}, & j' = 2j - 1 \\ 1, & j' = 2j \\ \frac{\Delta_k^{j'+1}}{\Delta_{k-1}^j}, & j' = 2j + 1 \\ 0, & \text{otherwise} \end{cases} \qquad (14)$$

$$g_k^{j'j} = \delta_{j',2j-1} . \qquad (15)$$

In one dimension, the element length from $x_k^{j-1}$ to $x_k^j$ is defined as $\Delta_k^j = x_k^{j'} - x_k^{j-1}$.

The action of $h_k^{j'j}$ applied at scale $k$ is to average, while $g_k^{j'j}$ acts to inject at scale $k$. This can be seen in Fig. 1 where the relationship between basis elements is given by

$$\phi_{k-1}^j = \sum_{j'} h_k^{j'j} \phi_k^{j'} \qquad (16)$$

$$\psi_k^j = \sum_{j'} g_k^{j'j} \phi_k^{j'} . \qquad (17)$$

With the Schauder basis defined, attention is turned to the question of stability. To begin the discussion on stability, several definitions are required.

**Definition 1** Given a basis $\mathbf{\Phi} = (\phi^1, \phi^2, \phi^3, \dots)^T$ for a subspace, $\mathcal{V}$, of a Hilbert Space, $\mathcal{H}$, with norm $\| \cdot \|_{\mathcal{H}} = \sqrt{(\cdot, \cdot)_{\mathcal{H}}}$, define the associated discretized matrix as

$$\mathbf{A}^{\mathbf{\Phi}} = a(\mathbf{\Phi}, \mathbf{\Phi})_{\mathcal{H}} . \qquad (18)$$

**Definition 2** Suppose the vector $\mathbf{\Phi}$ is a normalized basis for a Hilbert space, $\mathcal{H}$. Then $\mathbf{\Phi}$ is said to be stable in $\mathcal{H}$ if there exists constants $\alpha, \beta > 0$ such that

$$\alpha \|u\|_{\ell^2}^2 \le \|u^{\mathrm{T}} \mathbf{\Phi}\|_{\mathcal{H}}^2 \le \beta \|u\|_{\ell^2}^2, \qquad (19)$$

where $\|u\|_{\ell^2}^2 = \sum_i u_i^2$.

This definition of stability is equivalent to the definition of a Riesz basis with respect to the $\| \cdot \|_{\mathcal{H}}$ norm.

Now consider the following set of nested subspaces of $\mathcal{H}$,

$$\mathcal{V}_0 \subset \mathcal{V}_1 \subset \cdots \subset \mathcal{V}_k \cdots \subset \mathcal{H}, \qquad (20)$$

such that $\mathbf{\Phi}_k$ is a basis for $\mathcal{V}_k$.

**Definition 3** Given a nested sequence of subspaces of a Hilbert space, $\mathcal{H}$, and their respective bases $\{\mathbf{\Phi}_k\}_{k \in \mathbb{Z}^+}$ (normalized in $\mathcal{H}$), then $\{\mathbf{\Phi}_k\}_{k \in \mathbb{Z}^+}$ is uniformly stable in $\mathcal{H}$ if there exists constants $\alpha, \beta > 0$ independent of $k$ such that

$$\alpha \|u_k\|_{\ell^2}^2 \le \|u_k^{\mathrm{T}} \mathbf{\Phi}_k\|_{\mathcal{H}}^2 \le \beta \|u_k\|_{\ell^2}^2 \quad k = 0, 1, \dots \qquad (21)$$

With these definitions, the relationship between the stability of a basis and the condition number of an oper-

ator discretized in the basis is outlined. In specific, it will be demonstrated that a uniformly stable set of bases yields uniformly bounded condition numbers for the discrete operators independent of the refinement level.

**Lemma 1** Let $\mathbf{\Phi}_k$ be a basis for a subspace, $\mathcal{V}_k$, of a Hilbert space $\mathcal{H}$ such that $\mathbf{A}_k^{\mathbf{\Phi}}$ is positive definite. Let $\alpha = \|\mathbf{A}_k^{\mathbf{\Phi}^{-1}}\|_{\ell^2}^{-1}$ and $\beta = \|\mathbf{A}_k^{\mathbf{\Phi}}\|_{\ell^2}$. Then

$$\alpha \|u_k\|_{\ell^2}^2 \le \|u_k^{\mathrm{T}} \mathbf{\Phi}_k\|_{\mathcal{H}}^2 \le \beta \|u_k\|_{\ell^2}^2, \qquad (22)$$

and

$$\mathrm{cond}(\mathbf{A}_k^{\mathbf{\Phi}}) = \frac{\beta}{\alpha} . \qquad (23)$$

**Proof** The condition number of a symmetric matrix is given by $\mathrm{cond}(\mathbf{A}_k^{\mathbf{\Phi}}) = \|\mathbf{A}_k^{\mathbf{\Phi}}\|_{\ell^2} \|\mathbf{A}_k^{\mathbf{\Phi}}\|_{\ell^2}^{-1}$.

For the stability condition, note that

$$\frac{\|u_k^{\mathrm{T}} \mathbf{\Phi}_k\|_{\mathcal{H}}^2}{\|u_k\|_{\ell^2}^2} = \frac{(u_k^{\mathrm{T}} \mathbf{\Phi}_k, u_k^{\mathrm{T}} \mathbf{\Phi}_k)_{\mathcal{H}}}{u_k^{\mathrm{T}} u_k} = \frac{u_k^{\mathrm{T}} \mathbf{A}_k^{\mathbf{\Phi}} u_k}{u_k^{\mathrm{T}} u_k} . \qquad (24)$$

Moreover, for symmetric matrices

$$\beta = \|\mathbf{A}_k^{\mathbf{\Phi}}\|_{\ell^2} = \rho(\mathbf{A}_k^{\mathbf{\Phi}}) = \max_{u_k} \frac{u_k^{\mathrm{T}} \mathbf{A}_k^{\mathbf{\Phi}} u_k}{u_k^{\mathrm{T}} u_k} \qquad (25)$$

which gives us the upper bound. For the lower bound, observe that for symmetric positive definite matrices

$$\alpha = \|\mathbf{A}_k^{\mathbf{\Phi}^{-1}}\|_{\ell^2}^{-1} = \max_{u_k} \frac{u_k^{\mathrm{T}} \mathbf{A}_k^{\mathbf{\Phi}^{-1}} u_k}{u_k^{\mathrm{T}} u_k} \qquad (26)$$

The proof of Lemma 1 states that the best stability bounds are provided by the smallest and largest eigenvalues of the discretized matrix – a well known result.

Now consider the stability of a linear finite element basis, shown in Fig. 1a for $k = 2$, where $\mathcal{H} = L^2(0, 2)$ and $\mathcal{V}_k = \mathrm{span}\{\phi_k^1, j = 1, \dots, 2^{k+1} - 1\}$, $k = 0, 1, \dots, (27)$

$$\mathbf{\Phi}_k = (\phi_k^1, \phi_k^2, \dots, \phi_k^{2^{k+1}-1})^{\mathrm{T}} . \qquad (28)$$

**Lemma 2** The sequence of bases $\mathbf{\Phi}_k$ for the linear finite element basis is uniformly stable in $L^2(0, 2)$.

**Proof** Let the discretized matrix $\mathbf{A}_k^{\mathbf{\Phi}}$ be a diagonally scaled finite element mass matrix. In 1-D this is a symmetric tridiagonal matrix consisting of 1 on the diagonal and 1/4 on the non-zero off diagonals. Note that for the ensuing proof, $\mathbf{\Phi}_k$ has been normalized in a manner that is equivalent to diagonal preconditioning applied to the finite element mass matrix. Using Gershgorin's theorem, the eigenvalues for $\mathbf{A}_k^{\mathbf{\Phi}}$ lie in the interval $[1/2, 3/2]$ giving

$$\frac{1}{2} \|u_k\|_{\ell^2}^2 \le \|u_k^{\mathrm{T}} \mathbf{\Phi}_k\|_{\ell^2}^2 \le \frac{3}{2} \|u_k\|_{\ell^2}^2, \qquad (29)$$

and $\mathrm{cond}(\mathbf{A}_k^{\mathbf{\Phi}}) \le 3$. Therefore, the sequence of bases, $\mathbf{\Phi}_k$, is uniformly stable in $L^2$.

The stability associated with the finite element basis in $L^2$ is consistent with the empirical observation that the consistent mass matrix is well behaved in terms of its
condition number. In practice, this is reflected in the
ability to easily solve mass matrix dominated problems
with simple iterative techniques.

**Theorem 1** (Poincaré-Friedrich) For bounded domains, the $\mathcal{H}^m$ norm and semi-norm are equivalent in the sense that there exists constants $\alpha, \beta > 0$ such that

$$
\alpha|u|_{\mathcal{H}^m} \le \|u\|_{\mathcal{H}^m} \le \beta|u|_{\mathcal{H}^m} \quad \text{for } u \in \mathcal{H}_0^m . \quad (30)
$$

Recall that the norm and semi-norm for the Sobolev
space $\mathcal{H}^1$ for one space dimension are

$$
\|u\|_{\mathcal{H}^1} = \|u'\|_{L^2} + \|u\|_{L^2}
$$

$$
|u|_{\mathcal{H}^1} = \|u'\|_{L^2} .
$$

**Lemma 3** After re-normalization with respect to the $\mathcal{H}^1$
semi-norm, the sequence of bases $\mathbf{\Phi}_k$ is not uniformly
stable in $\mathcal{H}_0^1(0,2)$.

Proof Using the equivalence of the norm and semi-norm,
a sequence of vectors $u_k$ is constructed for which uniform
bounds are not possible. If

$$
u_k = \underbrace{(1, 1, 1, \dots, 1)}_{2^{k+1}-1} \quad \text{for } k = 0, 1, 2, \dots, \qquad (31)
$$

then

$$
\|u_k\|_{\ell^2}^2 = 2^{k+1} - 1
$$

and

$$
|u_k^T \Phi_k|_{\mathcal{H}^1}^2 = 1 \text{ for } k = 0, 1, 2, \dots, \quad (32)
$$

and stability for the basis can not be demonstrated.
Therefore, using the equivalence of the $\mathcal{H}^1$ norm and
semi-norm, there are no uniform bounds for the sequence
of bases $\Phi_k$ in $\mathcal{H}_0^1$.

Next, the hierarchical Schauder Basis is considered
where $\mathbf{\Phi}_0$ corresponds to the linear finite element basis
described above. Let

$$
\mathcal{W}_k = \operatorname{span}\{\phi_k^j, j = 0, 2, 4, \dots, 2^{k+1} - 2\} \quad (33)
$$

$$
\mathcal{V}_k = \mathcal{V}_0 \bigoplus_{j=1}^{k} \mathcal{W}_j \qquad (34)
$$

$$
\mathbf{\Psi}_k = (\phi_0^0, \phi_1^0, \phi_1^2, \dots, \phi_k^0, \phi_k^2, \dots, \phi_k^{2^{k+1}-2})^T. \quad (35)
$$

The hierarchical basis for $k = 2$ is shown in Fig. 1.

**Lemma 4** After re-normalization with respect to the $\mathcal{H}^1$
semi-norm, the sequence of bases $\mathbf{\Psi}_k$ for the Schauder basis
is uniformly stable in $\mathcal{H}_0^1(0,2)$.

Proof The $\mathcal{H}^1$ semi-norm of $\mathbf{\Psi}_k$ is equivalent to the $L^2$ norm of the orthogonal Haar wavelet basis. The Haar basis is an orthogonal basis for $L^2[0, 2]$ and thus $A_k^\Psi$ is the identity with respect to the semi-norm for $k = 0, 1, 2, \dots, \infty$.

**Lemma 5** After re-normalization with respect to the $L^2$
norm, the sequence of bases $\mathbf{\Psi}_k$ for the Schauder basis is
not uniformly stable in $L^2(0, 2)$.

*Proof* Define the sequence $u_k$ which gives the hat function centered at one with support $[1 - 2^k, 1 + 2^k]$ as

$$
u_k = (1, -1/2, -1/2, 0, -1/4, -1/4, \\
0, 0, 0, 0, -1/8, -1/8, 0, 0, 0, \dots) . \qquad (36)
$$

Thus,

$$
\|u_k\|_{\ell^2}^2 \ge 1, \quad (37)
$$

and

$$
\| u_k^T \mathbf{\Psi}_k \|_{L^2}^2 = \frac{2}{3} \frac{1}{2^k} . \qquad (38)
$$

Because $\|u_k^T \mathbf{\Psi}_k\|_{L^2}^2 \to 0$ as $k \to \infty$ the stability condition in Eq. (3) can not be met.

3

**The multi-scale finite element**

This section outlines the 1-D and 2-D multi-scale finite
elements. The multi-scale elements are based upon the
hierarchical Schauder basis with a change-of-basis incor-
porated at the element-level in order to make use of the
well-known finite element assembly procedure [6].

Before embarking on a description of the 1-D multi-
scale finite element, a brief review and interpretation of the
multi-scale transformation is presented. Recall from Eq.
(5) that $\mathbf{\Psi}_k = \mathbf{W}_k^T \mathbf{\Phi}_k$ where $\mathbf{W}_k$ is the multi-scale trans-
formation matrix. In order to make this transformation
concrete, consider the following example. Beginning with a
1-D grid consisting of 5 grid points and 4 linear finite
elements, the nodal basis will be decomposed into a
coarse-grid consisting of two elements and the associated
"psuedo-wavelets". This decomposition is shown sche-
matically in Fig. 2.

**Remark 1** The term pseudo-wavelet is used here to indi-
cate that the elements of $\mathcal{W}_k$ in the hierarchical Schauder
basis do not possess the property that their zeroth moment
is zero. However, the elements of $\mathcal{W}_k$ used to complete
the subspace $\mathcal{V}_{k-1}$ are uniformly stable in $\mathcal{H}^1$. Thus, the
term pseudo-wavelet seems appropriate.

<figure><img src="image_4.png" /><figcaption>Fig. 2. One-dimensional two-scale decomposition of a finite element nodal basis at scale-1</figcaption></figure>
In this example, the wavelet transform, $\Psi_1 = W_1^T \Phi_1$, is

$$ \left\{ \begin{matrix} \phi_0^3 \\ \psi_1^2 \\ \psi_1^4 \end{matrix} \right\} = \left[ \begin{matrix} 1/2 & 1 & 1/2 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{matrix} \right] \left\{ \begin{matrix} \phi_1^2 \\ \phi_1^3 \\ \phi_1^4 \end{matrix} \right\}, \quad (39) $$

where the subscript indicates the scale, and the superscript indicates the node number. From this, it is clear that the wavelet transform performs an averaging procedure to obtain the coarse-grid basis elements, $\Phi_0$, and an injection to obtain the coarse-grid pseudo-wavelets, $\Psi_1$.

The decomposition of nodal variables, $u$, may be accomplished with the wavelet transform, but its inverse is required to obtain the coarse-grid coefficients, i.e., $\Delta u = W^{-1}u$, where $\Delta u$ is the multi-scale component of the field. Relying on the inverse wavelet transform is impractical because the orthogonality constraint between the wavelets and scaling functions has been relaxed in the hierarchical Schauder basis. However, incorporating the hierarchical Schauder basis at the element level yields solution algorithms that compute the multi-scale solution directly and rely only on the reconstruction algorithm, i.e., $u = W\Delta u$. Thus, given a multi-scale representation of the field that corresponds to the multi-scale basis elements in Fig. 2, the reconstruction algorithm is simply

$$ \left\{ \begin{matrix} u_1^2 \\ u_1^3 \\ u_1^4 \end{matrix} \right\} = \left[ \begin{matrix} 1/2 & 1 & 0 \\ 1 & 0 & 0 \\ 1/2 & 0 & 1 \end{matrix} \right] \left\{ \begin{matrix} u_0^3 \\ \Delta u_1^2 \\ \Delta u_1^4 \end{matrix} \right\}. \quad (40) $$

As an aside, the wavelet transform is comprised of two components and both may be viewed in terms of a discrete convolution. Using the nomenclature introduced above for the wavelet transform, $T_1 = [H_1|G_1]$, where

$$ H_1 = \left\{ \begin{matrix} 1/2 \\ 1 \\ 1/2 \end{matrix} \right\}, \quad G_1 = \left[ \begin{matrix} 1 & 0 \\ 0 & 0 \\ 0 & 1 \end{matrix} \right]. \quad (41) $$

In this simple example of a two-scale decomposition, $W_1 = T_1$, but in general, the wavelet transform is computed recursively according to Eq. (10).

### The 1-D multi-scale element

The description of the multi-scale element begins with the linear finite element for which the shape functions are

$$ N_i = \frac{1}{2}(1 - \xi_i \xi). \quad (42) $$

Here, $\xi$ is the natural coordinate, $\xi_i$ is the nodal value of the natural coordinate, $i = 1, 2$, and $-1 \le \xi \le 1$.

The concept of scale is introduced at the element level by injecting degrees-of-freedom (DOF) that are supported by the pseudo-wavelets of the hierarchical Schauder basis. The 1-D multi-scale element is shown in Fig. 3 where a single “internal” degree-of-freedom located at $\xi = 0$ in the element is introduced at Scale-1, two DOF are introduced at Scale-2, and four DOF at Scale-3.

At Scale-1, the pseudo-wavelet is

$$ \psi_1(\xi) = \left\{ \begin{matrix} 1 + \xi & \text{if } -1 \le \xi \le 0 \\ 1 - \xi & \text{if } 0 \le \xi \le 1 \end{matrix} \right. . \quad (43) $$

More generally, the pseudo-wavelets for the multi-scale DOF may be written in terms of the translates and dilates of $\psi(\xi)$ as

$$ \psi_k(\xi) = \psi(\tilde{\xi}), \quad (44) $$

where

$$ \tilde{\xi} = 2^{k-1}(1 + \xi) - 2j - 1, \quad (45) $$

and

$$ \begin{aligned} 2^{2-k}j - 1 &\le \xi \le 2^{2-k}j + 2^{1-k} - 1 \\ 2^{2-k}j + 2^{1-k} - 1 &\le \xi \le 2^{2-k}j + 2^{2-k} - 1. \end{aligned} \quad (46) $$

Here $k$ indicates the scale, and $j$ indicates the translates in the element parametric space $(-1 \le \xi \le 1)$.

The derivatives of the shape functions yield constant functions that are orthogonal to the derivatives of the pseudo-wavelets at all scales. The derivatives of the pseudo-wavelets yield Haar wavelets as illustrated in Fig. 3, and at any given scale they are orthogonal with the derivatives at all other scales in the 1-D multi-scale element.

The reconstruction algorithm may be viewed as an element-by-element procedure that relies only on the multi-scale information in each element. The reconstruction is shown schematically in Fig. 3 where the DOF located at $\xi = -3/4$ is computed as a linear combination of the detail, $\Delta u$, at scales 0–3. Here, a canonical node-numbering scheme is used where the node numbers are $n_1, n_2$ at Scale-0, $n_3$ at Scale-1, $n_4, n_5$ at Scale-2, etc. as shown in Fig. 3.

With this numbering scheme, the reconstruction for the DOF located at $n_6$ may be written as

$$ u_6 = N_1 u_0^1 + N_2 u_0^2 + \psi_1^3 \Delta u_1^3 + \psi_2^4 \Delta u_2^4 + \psi_3^6 \Delta u_3^6, \quad (47) $$

where each of the basis elements is evaluated at $\xi = -3/4$. Thus, the reconstruction algorithm begins with the interpolant of the coarse-grid solution and injects refinements, or detail, up to the desired scale. The reconstruction may be written more generally as

$$ u_k^{2+2^{k-1}+j} = \sum_{i=1}^{N_{\text{npe}}} N_i u_0^i + \sum_{k=1}^{N_{\text{scale}}} \psi_k^{2+2^{k-1}+j} \Delta u_k^{2+2^{k-1}+j} \quad (48) $$

where the basis elements ($N_i$, and $\psi_k$) are evaluated at locations corresponding to the DOF location in the parent element.

### The 1-D multi-scale operators

The computation of the stiffness for the model problem, Eq. (1), is a straight-forward procedure that begins with the coarse-grid stiffness.

At the element level, i.e., Scale-0, this is simply

$$ K_0^e = \frac{\varepsilon}{h} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}, \quad (49) $$

where $h$ is the node-spacing for the coarse-grid. Making use of orthogonality, the diagonal stiffness entries for the scale DOF, i.e., for $k > 0$, are
<figure><img src="image_5.png" /><figcaption>Fig. 3. Basis elements and their piecewise derivatives for the one-dimensional multi-scale element with three refinement scales. The derivatives have been scaled by 1/k for scales k = 1, 2, 3</figcaption></figure>

$$K_k^{2+2^{k-1}+j} = \frac{2^{k+1}\epsilon}{h}. \qquad (50)$$

The resulting multi-scale element-level stiffness matrix is

$$K^e = \frac{\epsilon}{h} \begin{bmatrix} 1 & -1 & & & & & \\ -1 & 1 & & & & & \\ & & 4 & & & & \\ & & & 8 & & & \\ & & & & 8 & & \\ & & & & & \ddots & \\ & & & & & & 2^{k+1} \end{bmatrix}. \quad (51)$$

With this form of the element stiffness, only the coarse-scale terms contribute to the element assembly procedure since all internal DOF simply require a point evaluation and do not rely on information outside the element. Here, a uniform mesh refinement was applied to inject the scale DOF, but this is not a necessity for the multi-scale element.

The mass matrix computation is somewhat more involved because it requires the calculation of inner-products that involve the basis elements across scale. The computation of the multi-scale mass matrix consists of the following

$$ \begin{align} M_0^e &= \int_{\Omega^e} \gamma N_i N_j \, d\Omega, & M_{01} &= \int_{\Omega^e} \gamma N_i \psi_1 \, d\Omega, \nonumber \\ M_{10} &= \int_{\Omega^e} \gamma \psi_1 N_j \, d\Omega, & M_{11} &= \int_{\Omega^e} \gamma \psi_1 \psi_1 \, d\Omega \tag{52} \\ M_{12} &= \int_{\Omega^e} \gamma \psi_1 \psi_2 \, d\Omega, \dots \nonumber \end{align} $$

This series of element-level integrals leads to the finger-diagonal matrix structure shown in Fig. 4 with the inset coarse-grid element mass matrix.

A row-column lumping procedure is shown schematically in Fig. 4 for the multi-scale DOF at $n_7$ in the element. Note that the traditional row-sum lumping for the coarse-grid mass matrix can be used for the entries corresponding to the coarse-grid, but this type of mass lumping can not be used for the multi-scale DOF in the element. This concept is discussed further in the comparison between the multi-scale and linear finite elements.

### A multi-scale algorithm

As an example of a multi-scale algorithm, consider the model problem, Eq. (2), with $\epsilon = 1$, $\gamma = 0$, $f(x) = 1$, and essential boundary conditions, $u(0) = 0$, $u(1) = 0$. In
<figure><img src="image_6.png" /><figcaption>Fig. 4. Finger diagonal structure of multi-scale element mass matrix</figcaption></figure>

this example, three scale solutions are computed. Scale-0 corresponds to the coarse-grid solution using two elements. Scale-1 corresponds to the injection of one multi-scale DOF per element, while Scale-2 corresponds to the injection of two multi-scale DOF per element. This may be seen in Fig. 5 where the multi-scale DOF are shown relative to the elements ($e_1, e_2$) of the one-dimensional grid. At the finest resolution, the multi-scale solution corresponds to the solution associated with 8 linear finite elements. The multi-scale algorithm proceeds as follows.

**Algorithm 1 Multi-scale solution algorithm**

1. Form the coarse-grid operator $K_0$, and solve the coarse-grid problem, $K_0 U_0 = F_0$.

2. For each element, inject one scale DOF and solve for the wavelet coefficient, $\Delta u$.

$$ \{\Delta u_k\} = [K_k]^{-1}\{F_k\} \quad (53) $$

3. Compute the termination measure for the scale DOF injection. One possibility for the termination measure

<figure><img src="image_7.png" /><figcaption>Fig. 5. Exact solution and scale solutions for $k = 0, 1, 2$ (Node numbers 4–8 correspond to the insertion of multi-scale DOF)</figcaption></figure>

relies on stopping when the scale DOF are small relative to the overall solution. (Other stopping criteria have not been investigated.)

$$ \varepsilon = \frac{\Delta u_k^{2+2^{k-1}+j}}{\|u_k\|} \quad (54) $$

4. Repeat 1–3 until $\varepsilon$ is smaller than some user-specified criteria.

5. Perform the element-by-element multi-scale reconstruction using Eq. (48).

For the example problem, the multi-scale solution for $k = 1$ corresponds to $\Delta u_1^4 = \Delta u_1^5 = 1/32$, and for $k = 2$, the $\Delta u_2^6 = \Delta u_2^7 = 1/128$, and $\Delta u_2^8 = \Delta u_2^9 = 1/128$. After the reconstruction procedure, the scale DOF yield solution values that interpolate the exact solution – a result that is expected for linear problems. Similar results have been obtained for problems with non-linear source terms, and for problems with inhomogeneous essential and natural boundary conditions.

The algorithm presented for the 1-D multi-scale element possesses the property that all scale injection relies only on element-local data and does not require a re-solve of the coarse-grid problem to improve the solution. For problems that involve a significant mass-matrix contribution, the correction procedure outlined in the following sections is required when the row-column lumped mass is used.

**The 2-D multi-scale element**

Attention is now turned to the 2-D multi-scale element. As in the 1-D case, the 2-D bilinear element provides the element-level components of the global basis functions. To begin, Fig. 6 shows a four-patch of bilinear finite elements with the injected multi-scale DOF corresponding to $k = 1$. The configuration of the multi-scale DOF in the parent element is shown in Fig. 7. Like the shape functions, the components of the pseudo-wavelets at the element level take on a value of 1 at the DOF location, and they are zero at all other node locations.

<figure><img src="image_8.png" /><figcaption>Fig. 6. Four-patch of bilinear elements with multi-scale DOF</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 7a-f. Parent element and multi-scale DOF for $k = 1$. a reference element and basis elements for multi-scale DOF with $k = 1$, b $\psi_1^5$, c $\psi_1^6$, d $\psi_1^7$, e $\psi_1^8$, and f $\psi_1^9$</figcaption></figure>

The shape functions for the 2-D bilinear element are

$$N_i = \frac{1}{4}(1 + \xi_i \xi)(1 + \eta_i \eta), \quad (55)$$

where $i = 1, 2, 3, 4$, and $-1 \le \xi, \eta \le 1$. At the first scale, $k = 1$, the pseudo-wavelets are

$$\begin{align}
\psi_1^5 &= (1 - |\xi|)(1 - |\eta|), & \psi_1^6 &= \frac{1}{2}(1 + \xi)(1 - |\eta|), \nonumber \\
\psi_1^7 &= \frac{1}{2}(1 - |\xi|)(1 + \eta), & \psi_1^8 &= \frac{1}{2}(1 - \xi)(1 - |\eta|), \tag{56} \\
\psi_1^9 &= \frac{1}{2}(1 - |\xi|)(1 - \eta). & & \nonumber
\end{align}$$

In a more general way, the pseudo-wavelets for the multi-scale DOF may be written in terms of the translates and dilates of the basis functions at scale $k = 1$. The pseudo-wavelets in two-dimensions are

$$\psi_k^m(\xi, \eta) = \psi_1^m(\tilde{\xi}, \tilde{\eta}), \quad (57)$$

where

$$\begin{align}
\tilde{\xi} &= 2^{k-1}(1 + \xi) - 2j - 1, \tag{58} \\
\tilde{\eta} &= 2^{k-1}(1 + \eta) - 2j - 1, & & \nonumber
\end{align}$$
and

$$ 
\begin{aligned}
0 \le j \le k-1 \\
2^{2-k}j-1 \le \xi \le 2^{2-k}j+2^{2-k}-1 \\
2^{2-k}j-1 \le \tilde{\eta} \le 2^{2-k}j+2^{2-k}-1.
\end{aligned} 
\quad (59) 
$$

Here $k$ indicates the scale, $j$ indicates the translates in the element parametric space $(-1 \le \xi, \eta \le 1)$, and $m = 5, 6, 7, 8, 9$ for the element-local numbering of the pseudo-wavelets. With the basis elements defined this way, the use of recursion at the element level as indicated in Fig. 6 can be used to automate the computation of the mass and stiffness operators with a given scale of resolution, Nscale.

There are several key points regarding the 2-D multi-scale element. First, the orthogonality of the derivatives of the pseudo-wavelets is not preserved in two dimensions – even on an orthogonal grid. The finger-diagonal matrices that arise from this discretization can lead to excessive storage costs if the matrices are used without thresholding, row-column lumping, or element-by-element procedures. Surprisingly, the row-column lumping can be applied to both the multi-scale mass and stiffness operators, albeit only for the rows and columns corresponding to the multi-scale DOF. The 2-D multi-scale element is compatible with many h-adaptivity strategies being implemented in finite element codes today, however, a detailed discussion of this aspect is beyond the scope of this paper. The ability to use this element as a change-of-basis preconditioner is just beginning to be explored. Finally, the numerical performance of the 1-D and 2-D multi-scale element in terms of dispersive behavior, rate of convergence, etc. is identical to the linear (bilinear) element since any multi-scale solution

can be cast in terms of the reconstructed solution in the finite element basis at the finest grid scale.

## 4

### Multi-scale vs. linear elements

In this section, a comparison of the storage, matrix conditioning, and computational complexity for the hierarchical Schauder and linear finite element bases is presented in the context of the model elliptic problem developed in Sect. 2. The approximations inherent in using a lumped multi-scale mass matrix are discussed, and an efficient row-column lumping procedure with a lumped-mass correction algorithm is presented to illustrate that only the largest wavelet coefficients need be corrected in the approximate row-column lumped-mass solution. The row-column lumping procedure is applied to both the multi-scale mass and stiffness matrices in 2-D, and the lumped-mass correction algorithm is extended to account for the approximations introduced by the row-column lumped mass and stiffness matrices. Although a uniform discretization is considered in the ensuing discussion, it should be noted that in 1-D the hierarchical Schauder basis retains all of its properties for a non-uniform grid. However, in two-dimensions, orthogonality in the $H^1$ sense is lost even for uniform grids.

### 1-D comparison

The non-zero fill pattern for the mass and stiffness matrices are shown in Fig. 8 for both the hierarchical Schauder and the linear finite element bases. The formulae for the number of non-zeros associated with the mass and stiffness operators for both the hierarchical Schauder and finite element bases are shown in Table 1. Because the hierarchical Schauder basis diagonalizes the stiffness,

<figure><img src="image_10.png" /><figcaption>Fig. 8. The non-zero entries in the mass and stiffness matrices for the hierarchical Schauder and linear finite element bases for a mesh with 64 nodes</figcaption></figure>
<table><caption>**Table 1.** Formulae for the number of non-zeros in the 1-D mass and stiffness matrices where Ndof is the number of degrees-of-freedom, Ndof = $2^{k+1} - 1$, and k indicates the scale</caption><thead><tr><th>Basis - matrix</th><th>No. of non zeros</th></tr></thead><tbody><tr><td>Linear finite element - Mass matrix (M<sub>k</sub>)</td><td>3Ndof-2</td></tr><tr><td>Linear finite element - Stiffness matrix (K<sub>k</sub>)</td><td>3Ndof-2</td></tr><tr><td>Schauder - Mass matrix (M<sub>k</sub>)</td><td>(2k-1)<br>(Ndof + 1) + 3</td></tr><tr><td>Schauder - Stiffness matrix (K<sub>k</sub>)</td><td>Ndof</td></tr></tbody></table>

only Ndof (number of DOF) words of storage are required for the diagonal, although in practice a coarse-grid finite element matrix is required for the multi-scale element. The consistent mass matrix, although sparse, requires increased storage relative to the tri-diagonal finite element mass matrix due to the finger-diagonal structure. In comparison to the linear finite element, the storage for the multi-scale mass matrix scales with both the number of DOF, Ndof, and the scale, k.

As demonstrated in Sect. 2, the finite element mass and the hierarchical Schauder stiffness matrices are both well conditioned. In fact, the hierarchical Schauder stiffness is ideal in 1-D since the nodal change-of-basis results in a diagonal stiffness operator when the full wavelet transformation is applied as a preconditioner. In contrast, the finite element stiffness and the Schauder mass matrices are both poorly conditioned. To be more precise, the finite element basis is uniformly stable in $L^2$, but not in $\mathcal{H}^1$. However, the hierarchical Schauder basis is uniformly stable in $\mathcal{H}^1$, but not in $L^2$.

In order to illustrate the differences in stability between the finite element and the hierarchical Schauder bases, consider the condition numbers associated with the mass and stiffness matrices as shown in Table 2. Here, multiple levels of mesh refinement are considered with $k = 1$ corresponding to a mesh with 3 nodes (2 elements). The growth of the condition numbers for $K_k$ is seen to be proportional to $O(h^{-2})$ for the linear finite element bases, while the condition numbers for the mass matrix are bounded asymptotically at 3. The condition number associated with the combined finite element mass and stiffness ($M_k + K_k$) is dominated by the stiffness since $\varepsilon = \gamma = 1$ in this case. In contrast, the mass matrix for the hierarchical Schauder basis yields a condition number that grows approximately as $O(h^{-3/2})$ while the condition number for the stiffness is uniformly bounded at 1. The combined mass and stiffness operator ($M_k + K_k$) for the

Schauder basis is also uniformly bounded indicating the dominance of the stiffness.

Although the multi-scale mass matrix associated with the hierarchical Schauder basis is not as well conditioned as the finite element mass matrix, its primary drawback is the storage associated with the finger-diagonal fill pattern. However, the entries in the finger diagonals decrease in amplitude with increasing scale suggesting that the relative importance of the matrix entries become negligible with increasing scale. This effect is reflected in the integrals for the multi-scale mass matrix in Eq. (52). Therefore, either a mass lumping procedure or threshholding based on the relative size of the entries can be used to limit the matrix fill-in. In order to test this idea, an ad-hoc procedure for mass lumping is considered. In a finite element setting, mass lumping refers to the row-sum lumping procedure used to obtain a diagonal mass matrix. However, in the context of a multi-scale basis, the physical interpretation of mass lumping is not simple, and the idea is perhaps less well founded for the multi-scale element than for the linear finite element.

Regardless of these issues, experimentation with mass lumping strategies has suggested that a viable lumping procedure for the multi-scale mass matrix is to sum the values to the left and above the diagonal entry to the diagonal as shown in Fig. 4. The use of the row-column lumping procedure in Algorithm 1 permits the point-wise evaluation of the multi-scale DOF for problems where the mass matrix is required, but the subsequent reconstruction procedure can only yield an approximate solution. However, inspection of the resulting approximate solution, $\Delta u$, obtained with the lumped multi-scale mass indicates that a small correction to the largest wavelet coefficients can dramatically improve the solution quality. In other words, the row-column lumping procedure cannot be used over the entire range of scales without degrading the solution accuracy, i.e., the influence of mass lumping is “felt” more at the coarser grid scales. A lumped-mass correction algorithm has been developed to account for the approximate nature of the solution and proceeds as follows.

**Algorithm 2 Lumped-mass correction algorithm**

1. Row-column sum lump the multi-scale mass matrix resulting in a diagonal matrix, $\tilde{A} = M_k^L + K_k$. Here, the superscript L indicates a lumped operator. For the multi-scale element, only those entries in the mass matrix associated with the multi-scale DOF are lumped.

<table><caption>**Table 2.** Condition numbers for the 1-D mass and stiffness matrices after diagonal preconditioning for $\varepsilon = 1$, $\gamma = 1$. (Ndof = $2^{k+1} - 1$)</caption><thead><tr><th colspan="2">Mesh parameters</th><th colspan="3">Linear finite element</th><th colspan="3">Schauder</th></tr><tr><th>Scale (k)</th><th>Ndof</th><th>M<sub>k</sub></th><th>K<sub>k</sub></th><th>M<sub>k</sub> + K<sub>k</sub></th><th>M<sub>k</sub></th><th>K<sub>k</sub></th><th>M<sub>k</sub> + K<sub>k</sub></th></tr></thead><tbody><tr><td>1</td><td>3</td><td>2.09</td><td>6</td><td>4</td><td>7</td><td>1</td><td>1.23</td></tr><tr><td>2</td><td>7</td><td>2.71</td><td>25</td><td>18</td><td>18</td><td>1</td><td>1.27</td></tr><tr><td>3</td><td>15</td><td>2.92</td><td>103</td><td>74</td><td>49</td><td>1</td><td>1.29</td></tr><tr><td>4</td><td>31</td><td>2.98</td><td>414</td><td>295</td><td>121</td><td>1</td><td>1.29</td></tr><tr><td>5</td><td>63</td><td>3.00</td><td>1659</td><td>1181</td><td>288</td><td>1</td><td>1.29</td></tr><tr><td>6</td><td>127</td><td>3.00</td><td>6640</td><td>4724</td><td>667</td><td>1</td><td>1.29</td></tr><tr><td>7</td><td>255</td><td>3.00</td><td>26560</td><td>18900</td><td>1517</td><td>1</td><td>1.29</td></tr></tbody></table>
2. Calculate the approximate wavelet coefficients, $\Delta\tilde{u}$, by dividing the right-hand-side by the diagonal entries of $\tilde{A}$.

3. Permute the components of $\Delta\tilde{u}$ so that the N largest coefficients are at the top. Let $\Delta\tilde{u}^N$ be the largest N lumped-approximate wavelet coefficients and $\Delta\tilde{u}^R$ be the remaining coefficients such that

$$ \mathrm{perm}(\Delta\tilde{u}) = \begin{Bmatrix} \Delta\tilde{u}^N \\ \Delta\tilde{u}^R \end{Bmatrix}. \qquad (60) $$

4. Permute the rows and columns of $\tilde{A}$ to match the permutation of $\Delta\tilde{u}$ using a consistent multi-scale mass only for the $\Delta\tilde{u}^N$ unknowns. Here, the permuted matrix, $\tilde{A} = \tilde{M}_k + \tilde{K}_k$ is formed with the consistent mass matrix, $\tilde{M}_k$.

5. Solve the smaller system using the consistent mass matrix for $\Delta u^N$, i.e., the coarse-grid correction,

$$ A^N \Delta u^N = f^N - A^R \Delta \tilde{u}^R. \qquad (61) $$

Here, $A^N$ is the operator for the coarsest grid, and $A^R$ corresponds to the rows of $\tilde{A}$ associated with $\Delta\tilde{u}^R$ but formed with the consistent mass matrix. An efficient means of computing the matrix-vector product, $A^R\Delta\tilde{u}^R$ is through the use of the wavelet transform in Eq. (7). Multiplication of $\Delta u$ by the matrix $\tilde{W}$ is equivalent to the application of the wavelet transform, i.e., the reconstruction algorithm, which can be implemented in an $O(Nnp)$ operation algorithm. It should be noted that multiplication by $\tilde{W}^T$ is not equivalent to the inverse wavelet transform, i.e., $\tilde{W}^T\tilde{W} \neq I$. However, the properties of $\tilde{W}^T$ permit implementation of a matrix-vector product using convolution techniques similar to the reconstruction algorithm.

In practice, the DOF associated with the coarse-grid tend to be the largest and are the ones that need to be corrected. Thus, the permutation step can be replaced with a solve for the coarse-grid DOF followed by the injection of the multi-scale DOF. This algorithmic approach is exactly what the multi-scale finite element does since the coarse-grid solution is computed with only the injection of scale relying on the hierarchical Schauder basis. Figure 9 illustrates the effectiveness of the lumped-mass correction algorithm for the model problem with $\varepsilon = 1$, $\gamma = 1$, and $f = -x^8(x-2)^3$. In Fig. 9a, the exact solution is shown

with the approximate solution computed with the lumped multi-scale mass matrix. In Fig. 9b, the approximate solution has been improved by correcting only the three largest wavelet coefficients.

## 2-D comparison

In two dimensions, the behavior of the hierarchical Schauder basis changes rather significantly relative to the 1-D basis. The most apparent changes are the additional fill-in for the multi-scale mass matrix and the structure of the stiffness matrix which has a finger diagonal fill pattern rather than being diagonal. This is shown in Fig. 10 and is due to the multi-scale interaction of the tensor product basis in 2-D. The storage for the finger diagonal mass matrix grows as $4k^2N_{\mathrm{dof}}$ while the storage for the stiffness is proportional to $4kN_{\mathrm{dof}}$ where $N_{\mathrm{dof}} = (2^{k+1} - 1)^2$. In comparison, the storage for the finite element mass and stiffness matrices scales linearly in the number of DOF as $\sim 9N_{\mathrm{dof}}$. As in the 1-D case, the amplitude of the entries in the finger-diagonals of the multi-scale mass matrix decrease with increasing scale.

In addition to the change in the fill pattern, the $O(1)$ condition number for the 1-D Schauder stiffness becomes $O(\log(h^{-2}))$ in 2-D. Again, the condition number for the finite element stiffness grows as $O(h^{-2})$ regardless of the dimensionality. Table 3 illustrates how the condition number varies with mesh resolution for the model problem in Eq. (1). Here, the $L^2$ stability of the finite element basis is reflected in the uniformly bounded condition number for the mass matrix and the $O(h^{-2})$ growth of the condition number for the stiffness. In contrast, the condition number for the multi-scale mass matrix grows faster than $O(h^{-2})$ although firm estimates have not been made. For the stiffness, the condition number for $K_k$ grows as $\log(h^{-2})$. Although the condition number for the multi-scale mass matrix degrades rapidly, the condition number for $M_k + K_k$ is dominated by the stiffness and is much better behaved than the finite element basis as shown in Table 3.

Following the 1-D mass lumping, the concept of row-column lumping is applied in 2-D. However, here the lumping procedure is applied to both the multi-scale mass and the stiffness matrices. Although the row-sum lumping procedure is not valid for the finite element stiffness, the row-column lumping is possible for the multi-scale

<figure><img src="image_11.png" /><figcaption>Fig. 9a, b. Solutions, $u$, on $0 \le x \le 2$ obtained using the hierarchical Schauder basis with mass-lumping and a no corrections and b 3 wavelet coefficient corrections</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Table 3.** Condition numbers for the 2-D mass and stiffness matrices after diagonal scaling for $\varepsilon = 1$, $\gamma = 1$. (Ndof = $(2^{k+1}-1)^2$)</figcaption></figure>

<table><caption>**Fig. 10.** The non-zero entries in the mass and stiffness matrices associated with the hierarchical Schauder and linear finite element bases in 2-D for a 32 × 32 mesh</caption><thead><tr><th colspan="2">Mesh parameters</th><th colspan="3">Linear finite element</th><th colspan="3">Schauder basis</th></tr><tr><th>Scale (k)</th><th>Ndof</th><th>M<sub>k</sub></th><th>K<sub>k</sub></th><th>M<sub>k</sub> + K<sub>k</sub></th><th>M<sub>k</sub></th><th>K<sub>k</sub></th><th>M<sub>k</sub> + K<sub>k</sub></th></tr></thead><tbody><tr><td>1</td><td>9</td><td>4</td><td>3</td><td>3</td><td>49</td><td>7</td><td>8</td></tr><tr><td>2</td><td>49</td><td>7</td><td>13</td><td>11</td><td>334</td><td>18</td><td>19</td></tr><tr><td>3</td><td>225</td><td>9</td><td>52</td><td>43</td><td>2396</td><td>49</td><td>51</td></tr><tr><td>4</td><td>576</td><td>9</td><td>207</td><td>172</td><td>14621</td><td>121</td><td>123</td></tr><tr><td>5</td><td>1225</td><td>9</td><td>830</td><td>690</td><td>82832</td><td>288</td><td>290</td></tr></tbody></table>

representation of the stiffness – an artifact of the nodal change-of-basis. In other words, the multi-scale representation of the stiffness does not retain the “row-sum to zero” property of the nodal stiffness matrix. In order to illustrate this property, consider the application of Algorithm 2 in two-dimensions to the model problem with $\varepsilon = \gamma = 1$ and $f = xy(2-x)(2-y)$. However, in this case, the stiffness is replaced by the “lumped” stiffness, i.e., $\tilde{A} = M_k^L + K_k^L$.

Figure 11b and c show the lumped-approximate solution with no wavelet correction and the associated error relative to the exact solution. The approximate solution was obtained by performing one vector divide and a multi-scale reconstruction resulting in a very reasonable initial solution. In Fig. 11c, the peak in the error is associated with the largest approximate wavelet coefficient. Figure 11d illustrates the lumped-approximation after the largest wavelet coefficient has been corrected following Algorithm 2. Nearly an order of magnitude reduction in the error has been achieved by correcting only the largest wavelet coefficient. In Fig. 11e, the peaks in the error are associated with the next 8 largest wavelet coefficients. Although the row-column lumping procedure yields an approximate

solution, it is thought that this approach may generate a good initial solution for an iterative procedure or alternatively provide a multi-scale preconditioner.

Although the hierarchical Schauder basis is promising, the question of computational complexity remains. In 2-D, the change in fill pattern and matrix conditioning make computational complexity a more important issue in the comparison of the Schauder and linear finite element bases. In order to evaluate the computational complexity associated with the Schauder basis, a brief comparison of the cost of solving the model problem using a conjugate gradient method was performed. In this comparison, complete fill-in of the multi-scale mass and stiffness was accounted for ignoring the advantages of row-column lumping.

In order to perform the comparison, both the number of non-zero entries in the matrix, and the number of iterations required to solve the problem are required. Table 4 shows the number of non-zero entries and associated iteration count for the Schauder and linear finite element bases for the model problem with $0 \le \varepsilon \le 1000$ and $\gamma = 1$. For problems where the mass matrix dominates, $\varepsilon = 0$, the uniform stability of the finite element basis in $L^2$
<figure><img src="image_13.png" /><figcaption>a) Exact Solution: $u^e$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>b) Lumped Approximation ($u^h$) without Correction</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>c) Error: $u^e - u^h$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>d) Lumped Approximation ($u^h$) with 1 Correction</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>e) Error: $u^e - u^h$</figcaption></figure>

**Fig. 11a-e.** The exact and approximate solutions using the lumping procedure with and without wavelet correction for a

32 × 32 2-D mesh ($u^e$ is the interpolant of the exact solution, and $u^h$ is the discrete solution)

results in a computational cost that is bounded independent of the mesh resolution. In contrast, the hierarchical Schauder basis results in a lower-order computational cost for problems where the stiffness dominates due to the multi-scale preconditioning of the stiffness – despite the additional fill-in associated with the finger-diagonal matrices. This result is further emphasized by the data in

Table 5 for a “purely” elliptic problem with no mass matrix contribution.

In summary, Table 6 provides order of magnitude estimates for the computational complexity for the hierarchical Schauder and linear finite element bases. Here, the number of floating point operations required to solve the model problem grows as $9Ndof\sqrt{Ndof}$ for the finite
<table><caption>**Table 4.** Number of non-zeros and iteration count for the model problem using the **a** Schauder and **b** linear finite element bases in 2-D. The iteration count is based on using Jacobi preconditioned conjugate gradient with $0 \le \varepsilon \le 1000$, $\gamma = 1$</caption><thead><tr><th colspan="2">Mesh parameters</th><th rowspan="2">No. of non-zeros</th><th colspan="5">Iteration count</th></tr><tr><th>Scale (k)</th><th>Ndof</th><th>ε = 0</th><th>ε = 0.1</th><th>ε = 1</th><th>ε = 10</th><th>ε = 1000</th></tr></thead><tbody><tr><td colspan="9"><strong>a Schauder basis</strong></td></tr><tr><td>1</td><td>9</td><td>49</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td></tr><tr><td>2</td><td>49</td><td>729</td><td>11</td><td>10</td><td>10</td><td>9</td><td>9</td></tr><tr><td>3</td><td>225</td><td>6889</td><td>50</td><td>27</td><td>25</td><td>22</td><td>22</td></tr><tr><td>4</td><td>961</td><td>51529</td><td>187</td><td>48</td><td>42</td><td>37</td><td>37</td></tr><tr><td>5</td><td>3969</td><td>335241</td><td>488</td><td>73</td><td>61</td><td>52</td><td>51</td></tr><tr><td>6</td><td>16129</td><td>1990921</td><td>979</td><td>107</td><td>87</td><td>70</td><td>59</td></tr><tr><td colspan="9"><strong>b Linear finite element basis</strong></td></tr><tr><td>1</td><td>9</td><td>49</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td></tr><tr><td>2</td><td>49</td><td>361</td><td>6</td><td>6</td><td>7</td><td>7</td><td>7</td></tr><tr><td>3</td><td>225</td><td>1849</td><td>7</td><td>13</td><td>15</td><td>16</td><td>16</td></tr><tr><td>4</td><td>961</td><td>8281</td><td>6</td><td>27</td><td>31</td><td>31</td><td>31</td></tr><tr><td>5</td><td>3969</td><td>34969</td><td>5</td><td>54</td><td>62</td><td>63</td><td>63</td></tr><tr><td>6</td><td>16129</td><td>143641</td><td>3</td><td>110</td><td>125</td><td>126</td><td>126</td></tr></tbody></table>

<table><caption>**Table 5.** Storage requirements, matrix fill, and iteration count using Jacobi preconditioned conjugate for the Schauder basis with the “purely” elliptic model problem with $\varepsilon = 1$, $\gamma = 0$</caption><thead><tr><th colspan="2">Mesh parameters</th><th rowspan="2">No. of non-zeros</th><th rowspan="2">Iteration count</th></tr><tr><th>Scale (k)</th><th>Ndof</th></tr></thead><tbody><tr><td>1</td><td>9</td><td>33</td><td>3</td></tr><tr><td>2</td><td>49</td><td>329</td><td>9</td></tr><tr><td>3</td><td>225</td><td>2265</td><td>22</td></tr><tr><td>4</td><td>961</td><td>13113</td><td>37</td></tr><tr><td>5</td><td>3969</td><td>68985</td><td>51</td></tr><tr><td>6</td><td>16129</td><td>342265</td><td>59</td></tr></tbody></table>

<table><caption>**Table 6.** Computational complexity for the linear finite element and the hierarchical Schauder basis in 2-D where $N_{dof} = (2^{k+1}-1)^2$, and $k$ indicates the scale</caption><thead><tr><th>Basis</th><th>No. of non-zeros</th><th>Iteration count</th><th>Floating point operations</th></tr></thead><tbody><tr><td>Linear finite element (M<sub>k</sub> + K<sub>k</sub>)</td><td>9N<sub>dof</sub></td><td>√N<sub>dof</sub></td><td>9N<sub>dof</sub>√N<sub>dof</sub></td></tr><tr><td>Schauder basis (M<sub>k</sub> + K<sub>k</sub>)</td><td>4k<sup>2</sup>N<sub>dof</sub></td><td>12k</td><td>48k<sup>3</sup>N<sub>dof</sub></td></tr><tr><td>Schauder basis (K<sub>k</sub>)</td><td>4kN<sub>dof</sub></td><td>12k</td><td>48k<sup>2</sup>N<sub>dof</sub></td></tr></tbody></table>

element basis. In comparison, the number of floating point operations required to solve $M_k + K_k$ using the Schauder basis scales as $48k^3N_{dof}$ and as $48k^2N_{dof}$ for the purely elliptic problem. From this data, the Schauder basis does eventually have a lower computational cost, but the benefit is only apparent for relatively large problems. For the “purely” elliptic operator, $K_k$, the hierarchical Schauder basis wins, i.e., has lower computational cost, for mesh resolution exceeding $N_{dof} = 10^5$. However, for $M_k + K_k$, the Schauder basis does not win until the mesh resolution exceeds $N_{dof} \ge 10^8$.

# 5

## Summary and conclusions

In the search for a wavelet basis with good numerical and computational performance for multi-scale simulations, the following “shopping list” of characteristics was

developed as the goal for the ideal wavelet basis. The attributes include: (1) Low order, e.g., linear, for computational efficiency, (2) Consistent reproduction of polynomials, e.g., reproduce $\{1, x, y, xy\}$ in two-dimensions, (3) Nodal, i.e., possesses the Kronecker delta property, (4) Hierarchical: $\mathcal{V}_1 = \mathcal{V}_0 \oplus \mathcal{W}_1$, (5) Element based – compatible with isoparametric finite elements, (6) Analytic expressions for the basis elements $\phi$ and $\psi$, (7) Easy treatment of boundary conditions, (8) Good numerical performance, e.g., dispersion characteristics, truncation error, etc., (9) Appropriate for both Eulerian and Lagrangian computations, (10) Computationally efficient decomposition and reconstruction of fields, and (11) Extensible to multiple spatial dimensions. Based on these characteristics and the results of this exploratory effort the following conclusions are drawn.

The hierarchical Schauder basis, and the multi-scale elements, are prototypical of what the ideal multi-scale basis should be although this basis does not satisfy all 11 characteristics of an ideal basis. Unfortunately, the storage and computational cost associated with the finger-diagonal matrices from this type of basis is a significant penalty. However, the use of ad-hoc lumping procedures ameliorates this problem and offers the potential for the development of fast, simple preconditioners. Currently, the real value of the multi-scale elements lies in the application to elliptic problems. A preconditioner based on the multi-scale element would be of great value in applications where a dominant elliptic component is present, e.g., in time-dependent incompressible flow and quasi-static electro-magnetics.

As demonstrated in the discussion of the hierarchical Schauder basis, it is difficult to construct a wavelet basis that is stable in both $L^2$ and in $\mathcal{H}^1$, i.e., for all possible combinations of mass and stiffness operators. The design of wavelet bases that are customized for a specific partial differential equation remains an open topic of active research, and the use of wavelet bases for the solution of partial differential equations remains a research topic that is currently centered in the mathematics community. However, hierarchical solution procedures that use wavelets tailored to the physical problem appear to be the most viable candidates for using wavelet bases.
References

[1] Canuto C, Tobaco A, Urban K (1998) The wavelet element method part ii: Realization and additional features in 2d and 3d, Tech. Rep. 98-17, Institute for Computer Applications in Science, NASA Langley Research Center, Hampton, VA April 1998

[2] Christon MA, Baty RS, Burns SP, Roach DW, Trucano TG, Voth TE, Weatherby JR, Womble D (1998) An investigation of wavelet bases for grid-based multi-scale simulations – Final Report, Tech. Rep. SAND 98-2456, Sandia National Laboratories, Albuquerque, New Mexico, September 1998

[3] Dahmen W (1997) Wavelet and multiscale methods for operator equations. Acta Numerica 55–228

[4] Dahman W, Micchelli CA (1995) Biorthogonal wavelet expansions, Tech. Rep. 114, Institute for Geometry and Applied Mathematics – RWTH Aachen, Aachen, Germany, May 1995

[5] Grossmann A, Morlet J (1984) Decomposition of hardy functions into square integrable wavelets of constant shape. SIAM Journal on Mathematical Analysis 15:723–736

[6] Hughes TJR (1987) The Finite Element Method, Prentice-Hall, Inc., Englewood Cliffs, New Jersey

[7] Meyer Y (1993) Wavelets Algorithms and Applications, SI-AM, Philadelphia, Pennsylvania

[8] Morlet J (1983) NATO ASI Series, Issues in Acoustic Signal/ Image Processing and Recognition, vol. 1, Springer, Berlin, 1983, ch. Sampling theory and wave propagation, pp. 233–261

[9] Morlet J, Arens G, Fourgeau I, Giard D (1982) Wave propagation and sampling theory, Geophys. 47:203–236

[10] Yserentant H (1986) On the multi-level splitting of finite element spaces. Numerishce Mathematik 49:379–412

[11] (1990) Two preconditioners based on the multi-level splitting of finite element spaces. Numerische Mathematik 58:163–184

[12] Zienkiewicz OC (1994) The Finite Element Method, McGraw-Hill, England, 4th ed., 1994, ch. 7, pp. 110–149

[13] Zienkiewicz OC, Kelley DW, Gago J, Babuska I (1982) The mathematics of finite elements and applications IV, Academic Press, London, England, 1982, ch. Hierarchical finite element approaches, error estimates and adaptive refinement