Wavelet Based Numerical Homogenization

Bjorn Engquist

Abstract. In analytic homogenization, a differential equation and its solution with multiple scales are replaced by an approximating equation and its corresponding smoother solution with fewer scales. The scales related to the shortest wavelengths are eliminated. We shall start from a discretization of the original differential equation, which includes all the scales. The solution and the difference operator will be represented in a wavelet basis and the homogenized discrete operator will correspond to a particular form of an approximative projection onto the coarser scales. We shall show that this new operator inherits many of the properties of the original discrete operator, including sparseness. Some numerical examples will be presented and comparisons with the analytic homogenization process will be given. We shall also discuss direct coarse grid approximation.

1. Introduction. Homogenization is a classical analytical way to approximate the effect of some classes of periodic or stochastic oscillations. The problem is often formulated as follows. Consider a set of operators L_{ϵ}, indexed by the small parameter ϵ, and a right hand side f. Find the homogenized operator L̄ defined by

L_{ϵ}u_{ϵ} = f,  lim_{ϵ→0} u_{ϵ} = ū,  L̄ū = f.  (1)

In certain cases the convergence above and existence of the homogenized operator can be proved, [3].

In the d-dimensional elliptic case, let A(y) ∈ R^{d×d} be one-periodic in each of its arguments and let I_{d} denote the unit square. It can then be shown, [3], that

L_{ϵ} = −∇ · (A (\frac{x}{ϵ}) ∇) , L̄ = −∇ · (H∇), H = ∫_{I_{d}} A(y) − A(y)Dχ(y)dy, (2)

where Dχ is the Jacobian of the vector valued function χ(y) ∈ R^{d}, whose components χ_{k} are given by solving the so called cell problems

∇ · (A(y)∇χ_{k}) = ∑_{i=1}^{d} \frac{∂}{∂y_{i}} a_{ik}(y),  k = 1, …, d,  (3)

with periodic boundary conditions for χ_{k}. Note that H is a constant matrix. See [9] for a direct numerical application of this analytic formalism.

In this paper we present a general procedure for constructing numerical subgrid models to be used on a coarse grid where the smallest scales are not resolved. As in analytic homogenization the subgrid phenomina can be oscillations. The wave length ϵ in the oscillations may be smaller than the the typical grid step size h. The objective is to find models that accurately reproduce the effect of subgrid scales and that in some sense are similar to the original differential operator as is the case in analytic homogenization. The starting point is a finite-dimensional approximation, Lu = f, of a differential equation where L approximates the differential operator and u the solution. The operator L can be written on the form

L = P (Δ, A, h, ϵ),  (4)

where Δ is a collection of difference operators, A are discretized variable coefficients, typically diagonal matrices, h represents the grid size, and ϵ the smallest scale of significance in the problem.

We shall first briefly discuss the possibility of directly discretizing (4) on a coarse grid, h > ϵ. In general, for finite difference and finite element methods, a reasonable number of grid points or elements are required per wave length of the oscillation, h << ϵ. Phase and group velocity errors will otherwise be O(1).

For a special type of problems and numerical methods it is, however, possible to prove convergence in a weak sense even if the oscillations are not resolved on the computational grid. These types of techniques are studied in [10], [11] and commented on in section 2.

For the wavelet based homogenization technique we start with a resolved discretization, h << ϵ, and a coarse grid approximation. The specific scale ϵ does not play a role any longer and is dropped in the notation.

We seek a finite dimensional operator L̃ and a right hand side f̄ with the following properties. First, L̃ũ = f̄and ũ is a projection of u onto a lower dimensional subspace. Second, L̃ can be written on the same form as L,

L̃ = P (Δ, H, h̄),  (5)

but with h̄ >> h and the structure of H close inheting essential properties from the structure of A, typically diagonal dominance and sparsity. The sparsity of the discrete operator is important and corresponds to L̃ being an approximation of a differential operator. We interpret H as the subgrid model of A. If A corresponds to a material coefficient, H can be seen as the effective material coefficient. The procedure outlined above resembles that of the analytic homogenization technique used for the continuous case, see section 4. In view of this, we will call L̃ the homogenized operator. See Bensoussan et al., [3], for a thorough presentation of classical homogenization.

Our method is based on multiresolution analysis with wavelet projections and approximation of the discrete operator. Although it can be used with any type of discretization, it is algebraic and, in the present form, only deals with linear systems of equations. The great advantage of this procedure to derive subgrid models is its generality. It can be used on any system of differential equations and does not require separation into the distinct O(ϵ) and O(1) scales or periodic coefficients. It can also be used to test if it is physically reasonable to represent fine scale effects on a coarse grid with a local operator.

This work was initially presented in Dorobantu and Engquist [8], Andersson, Engquist, Ledfelt and Runborg [1], and based on ideas from Brewster and Beylkin, [5]. See also [13] for analysis in the one-dimensional case. Moreover, there are similarities with numerical homogenization based on techniques from algebraic multigrid, [15,16] and from the use of special purpose finite element methods, [14].

2. Direct discretization. Let us first consider the simple approach of using a coarse grid even if not all scales of the original differential equation are clearly resolved. For solutions which are highly oscillatory relative to the grid discretization, numerical techniques without phase velocity errors are needed. In [10], [11] particle scheme or method of characteristics approximations of hyperbolic partial differential equations are analyzed. For a restricted class of schemes it is possible to prove convergence, or weak convergence, in L^{p} of the numerical approximation to the analytic solution as h → 0 essentially independent of ϵ. Convergence essentially independent of ϵ means that a set of ratios of h/ϵ with arbitrary small Lebesque measure must be excluded to avoid resonance, [10], [11].

One simple but typical example for which a rigorous theory is possible is the method of characteristics for the Carleman equations,

\frac{∂u}{∂t} + \frac{∂u}{∂x} + u^{2} − v^{2} = 0

\frac{∂v}{∂t} − \frac{∂v}{∂x} + v^{2} − u^{2} = 0  (6a)

u(x, 0) = a(x, x/ϵ) v(x, 0) = b(x, x/ϵ)

a(x, y), b(x, y), 1-periodic in y
u(x_{j}, t_{n}) ∼ u_{j}^{n}
x_{j} = jΔx, t_{n} = nΔt, Δt = Δx,

u_{j}^{n+1} = u_{j−1}^{n} + Δt((v_{j−1}^{n})^{2} − (u_{j−1}^{n})^{2}),
v_{j}^{n+1} = v_{j+1}^{n} + Δt((u_{j+1}^{n})^{2} − (v_{j+1}^{n})^{2}),  (6b)
u_{j}^{0} = a(x_{j}, x_{j}/ϵ)
v_{j}^{0} = b(x_{j}, x_{j}/ϵ)

The homogenization theory of Tartar [17] applies to the differential equations (6a) and is also used in the convergence proof. The local truncation errors are large for h > ϵ and a cancelation of the errors must be established. The theorem gives strong convergence in L_{∞} essentially independent of ϵ as h → 0.

The wavelet based type of homogenization was derived in order to handle wider classes of differential equations.

3. Wavelet based homogenization. Given the full discrete solution operator on a fine grid we wish to find an operator of lower dimension that extracts only the coarse scales of the solution. Let V_{j} and W_{j} refer to the usual scaling and wavelet spaces, see e.g. [7]. Then, for a solution in V_{j+1} = V_{j} ⊕ W_{j}, the coarse scale is represented by V_{j}, and we are thus interested in the operator that yields the solution’s projection onto V_{j}.

Consider the equation

L_{j+1}U = F,  U, F ∈ V_{j+1},  (7)

originating from a discretization of a differential equation, where U, in the Haar case, is identified as a piecewise constant approximation. We introduce the orthogonal transformation

W_{j} : V_{j+1} → W_{j} × V_{j}, W_{j}U ≡ [ \begin{matrix}U_{h} \\ U_{l} \end{matrix} ] U_{h} ∈ W_{j}, j U_{l} ∈ V_{j},  (8)

and note that the linear operator W_{j}L_{j+1}W_{j}^{T} can be decomposed into four operators L_{j+1} = A_{j} + B_{j} + C_{j} + L_{j}, acting between the subspaces V_{j} and W_{j}, and such that (7) becomes

[ \begin{matrix}A_{j} & B_{j} \\ C_{j} & L_{j} \end{matrix} ] [ \begin{matrix}U_{h} \\ U_{l} \end{matrix} ] = [ \begin{matrix}F_{h} \\ F_{l} \end{matrix} ] , U_{h}, F_{h} ∈ W_{j}, U_{l}, F_{l} ∈ V_{j}.  (9)

when we apply W_{j} from the left. Block Gaussian elimination now gives an equation for U_{l}, the coarse part of the solution,

L̄_{j}U_{l} = F̄_{j}, L̄_{j} = L_{j} − C_{j}A_{j}^{−1}R_{j}, F̄_{j} = F_{l} − C_{j}A_{j}^{−1}F_{h}.  (10)

Hence, our new “coarse grid operator” L̄_{j} is the Schur complement of W_{j}L_{j+1}W_{j}^{T}. We also get the homogenized right hand side, F̄_{j}.

For higher dimensions, a standard tensor product extension of the multiresolution analysis allows us to use essentially the same derivation as above to obtain coarse grid operators.

We should note that in general L̄_{j} will not be sparse even if L_{j+1} is. For the method to be efficient we must be able to approximate L̄_{j} with a sparse matrix L̃_{j}. This is possible in many important cases. The fact that L̄_{j} is approximately sparse is fundamental. The finite dimensional operator L̃_{j} is our numerically homogenized operator.

The homogenization procedure can be applied recursively on L̄_{j} to get L̄_{j−1} and so on. This can easily be verified when L_{j+1} is symmetric positive definite. Furthermore, the condition number will not deteriorate. From [2], Chapter 3 with L_{j+1} = L_{j+1}^{T} and

c_{1}∥U∥^{2} ≤ ⟨L_{j+1}U, U⟩ ≤ c_{2}∥U∥^{2},  ∀U ∈ R^{2^{j+1}}  (11)

we have the same constants c_{1}, c_{2},

c_{1}∥V ∥^{2} ≤ ⟨L_{j}V, V ⟩ ≤ c_{2}∥V ∥^{2}, ∀V ∈ R^{2^{j}},  (12)

where L̄_{j} is defined by (10) and ⟨u, v⟩ = ∑_{k} ū_{k}v_{k}. For the first step in the process an improvement in the condition number can often be estimated from

⟨L̄_{j}V, V ⟩ = ⟨(L_{j} − B_{j}^{T} Λ_{j}^{−1}B_{j})V, V ⟩ = ⟨L_{j}V, V ⟩ − ⟨A_{j}^{−1}B_{j}V, B_{j}V ⟩
≤ ⟨L_{j}V, V ⟩.  (13)

When the operator L_{j+1} is derived from a finite difference, finite element or finite volume discretization, it is sparse and of a certain structure. In one dimension it might for instance be tridiagonal. However, as remarked above, the matrix L̄_{j} is not sparse since A_{j}^{−1} is usually dense. Computing all components of L̄_{j} would be inefficient. Fortunately, L̄_{j} will be diagonal dominant in many important cases. For instance, in [8] we proved that for a class of elliptic problems the matrix elements of L̄_{j} decay exponentially away from the diagonal. We are then able to find a sparse matrix that is a close approximation of L̄_{j}. In general the sparse approximation property follows from the analysis of Calderon-Zygmund operators in Beylkin, Coifman and Rokhlin, [4].

One simple way approximate L̄_{j} is to set all components outside a prescribed bandwidth equal to zero. Let us define truncation of M to bandwidth ν as

trunc(M, ν)_{ij} = \{ \begin{matrix}M_{ij}, & if 2|i − j| ≤ ν − 1 \\ 0, & otherwise. \end{matrix}  (14)

There are natural extensions to multi dimensions. This procedure was introduced in [4] and used in [8]. We propose that L̄_{j} be projected onto banded form in a more effective manner. Let {v_{j}}_{j=1}^{ν} be a set of linearly independent vectors in R^{2^{j}}. We define the band projection, band(M, ν), of a matrix M as the projection of M onto the subspace of matrices with bandwidth ν such that

Mx = band(M, ν)x,  ∀x ∈ span{v_{1}, v_{2}, …, v_{ν}}.  (15)

In our setting M will usually operate on vectors representing smooth functions, for instance solutions to elliptic equations, and a natural choice for the v_{j} vectors in one dimension are thus the first ν polynomials,

v_{j} = {1^{j−1}, 2^{j−1}, …, N^{j−1}}^{T} .  (16)

For the case ν = 1 we should remark that we get the standard “masslumping” of a matrix.

This technique is similar to the probing technique used by Chan et al, [6]. In that case the vectors v_{j} are sums of unit vectors. Other probing techniques have been suggested by Axelsson, Pohlman and Wittum, see e.g. Chapter 8 in [2]. In some cases the band projection technique only gives improvements for small values of ν, see figure 1. Numerical evidence indicate that for small values of ν, the band projection technique is quite efficient.

In figure 1 the results from different types of compressions are given for a simple test example. The differential equation is

\frac{d}{dx}a(x)\frac{d}{dx}u(x) = 1, 0 ≤ x ≤ 1,  (17)

u(0) = 0, \frac{d}{dx}u(1) = 0,

and it is approximated by centered differences. The variable coefficient a(x) > 0 is a uniformly distributed random function. It is clear that the divergence form of the operator with an explicit H-matrix is important. Reducing the number of non-zero elements in H is more efficient than in L̄_{j}.

Figure 1: The effect of different types of projections of L_{j} in the

homogenized approximations of (17).

For simple elliptic problems with periodic coefficients it is possible to prove rigorous error estimates for the projection techniques, [1], [8]. The theorems typically state that a finite approximation error δ can be achieved by an operator L̃_{j}, the matrix of which has only a low number (m) of non-zero elements in each row. The number m depends on δ but at most as log h on h.

4. Comparisons and examples. There is a striking relationship between the analytically homogenized operator in (2) and the Schur complement L̄_{j} in (10). The first term in both the expressions represent averaged operators, L_{j} in a discrete sense and ∫_{I_{d}} A is an integral sense. In both formulations a correction term for the high frequency interaction is subtracted from the average. Furthermore, in the correction term χ is the solution of an elliptic equation and A_{j}^{−1} is an analogous discrete positive definite solution operator.

In the classical analytical setting (2) homogenization gives the asymptotic expansion

u_{c}(x) = ū(x) + ϵu_{1}(x, x/ϵ) + O(ϵ^{2})  (18)

for the solution, see [3]. The techniques described in [14], [15] and [16] give numerically homogenized operators with coarse grid solutions which directly samples (18). The oscillatory term ϵu_{1}(x, x/ϵ) is contained in the solution. In the wavelet homogenization the solution is a projection of u_{ϵ} onto a coarse scale V_{j} space. The influence of the u_{1} term, which oscillates with mean zero [4], is then significantly reduced in orders of ϵ, depending on the type of wavelet basis. Another advantage of using wavelets is the favorable compression ratio for A_{j}^{−1}, [4].

Figure 2: coefficient a(x) in the Helmholtz equation

A number of problems have been studied numerically. Some are simple test cases in one space dimension for the analysis of different properties of the wavelet based homogenization technique. One example was given in section 3 above. See also [1], [7] and [8]. Others originate from more realistic simulations of e.g. wave propagation in fiber optics and flow in porous media, [12].

A practical application is subcell modelling in the form of the coarsening of a mesh refinement. Sometimes mesh refinements are necessary in order to resolve small geometric details. The refined mesh increases the computational complexity. The numerical homogenization can produce an operator on a uniform coarser grid which inherits the correct resolution from the finer grid. In figure 2 and 3 an example of such a grid coarsening is given.

A standard centered finite difference approximation of the Helmholtz equation,

−∇a(x)∇u − k^{2}u = 0, u : R^{2} → R,  (19)

is projected onto a coarser wavelet space. The figure displays different levels of homogenizations, and it is clear that the quality of the

The application comes from the study of wave propagation through a slit for the analysis of electromagnetic compatibility. The figure 4 gives the structure of the matrix corresponding to L̃_{j}.

Figure 3: Discrete solution of the Helmholtz equation for different

levels of numerical homogenization

Figure 4: Non-zero elements in the matrix representation of the projected

L̄_{j} corresponding to the solution given in figure 3.

References

[1] U. Andersson, B. Engquist, G. Ledfelt and O. Runborg. A contribution to wavelet-based subgrid modeling. To appear.

[2] O. Axelsson. Iterative Solution Methods. Cambridge University Press, 1991.

[3] A. Bensoussan, J.-L. Lions, and G. Papanicolau. Asymptotic Analysis for Periodic Structures. North-Holland Publ. Comp., The Netherlands, 1978.

[4] G. Beylkin, R. Coifman, and V. Rokhlin. Fast wavelet transforms and numerical algorithms I. Comm. Pure Appl. Math., 44:141-183, 1991.

[5] M. Brewster and G. Beylkin. A multiresolution strategy for numerical homogenization. Applied and Computational Harmonic Analysis, 2:327-349, 1995.

[6] T. Chan and T. Mathew. The interface probing technique in domain decomposition. SIAM J. Matrix Anal. Appl., 13(1):212-238, January 1992.

[7] I. Daubechies. Ten Lectures on Wavelets. SIAM, 1991.

[8] M. Dorobantu and B. Engquist. Wavelet-based numerical homogenization. SIAM J. on Num. Anal., 35(2):540-559, April 1998.

[9] L. Durlofsky. Numerical calculation of equivalent grid block permeability tensors for heterogenuous porous media. Water Resour. Res., 27:699-708, 1991.

[10] B. Engquist and T. Hou. Particle method approximation of oscillatory solutions to hyperbolic differential equations. SIAM J. Num. Analysis, 26:289-319, 1989.

[11] B. Engquist and J.-G. Liu. Numerical methods for oscillatory solutions to hyperbolic problems. Comm. Pure Appl. Math., 46:1-36, 1993.

[12] B. Engquist, P.-O. Persson and O. Runborg. Simulation of optical resonance filters using wavelet based homogenization. To appear.

[13] A.C. Gilbert. A comparison of multiresolution and classical one-dimensional homogenization schemes. Appl. Comput. Harmon. Anal., 5(1):1-35, 1998.

[14] T.Y. Hou and X.H. Wu. A multiscale finite element method for elliptic problems in composite materials and porous media. J. Comput. Phys., 134(1):169-189, 1997.

[15] S. Knapek. Matrix-dependent multigrid-homogenization for diffusion problems. SIAM J. Sci. Stat. Comp., 1998. To appear.

[16] N. Neuss. Homogenisierung und Mehrgitter. PhD thesis, Fakultät Mathematik Universität Heidelberg, 1995.

[17] L. Tartar, Etudes des oscillations dans les equations aux derivées partielle non linéaires. Lecture Notes in Physics, 195:384-412, Springer-Verlag, 1984.

Bjorn Engquist

Department of Mathematics, UCLA, Los Angeles, CA 90095-1555 USA and NADA, KTH, 10044 Stockholm, Sweden