Wavelet Based Fast Solution of Boundary Integral Equations

This paper presents a wavelet Galerkin scheme for the fast solution of boundary integral equations. Wavelet Galerkin schemes employ appropriate wavelet bases for the discretization of boundary integral operators which yields quasi-sparse system matrices. These matrices can be compressed such that the complexity for solving a boundary integral equation scales linearly with the number of unknowns without compromising the accuracy of the underlying Galerkin scheme. Based on the wavelet Galerkin scheme we present also an adaptive algorithm. By numerical experiments we provide results which demonstrate the performance of our algorithm. AMS Subject Classification: 47A20, 65F10, 65F50, 65N38, 65R20.


Introduction
Many mathematical models concerning e.g.field calculations, flow simulation, elasticity or visualization are based on operator equations with global operators, especially boundary integral operators.Traditional discretizations will then lead in general to possibly very large linear systems with densely populated matrices.Therefore, the complexity for solving such equations is at least O(N 2 J ), where N J denotes the number of equations.This fact restricts the maximal size of the linear equations seriously.
Modern methods for the fast solution of boundary integral equations reduce the complexity to a suboptimal rate, i.e.O(N J log α N J ), or even an optimal rate, i.e.O(N J ). Prominent examples for such methods are the fast multipole method [21] and the panel clustering [22].Wavelet compression [1] offers a further tool for the fast solution of integral equations.In fact, a Galerkin discretization with wavelet bases results in quasi-sparse matrices, i.e., the most matrix coefficients are negligible and can be treated as zero.Discarding these nonrelevant matrix coefficients is called matrix compression.It has been shown in [9,31] that only O(N J ) significant matrix coefficients remain.
Concerning boundary integral equations, a strong effort has been spent on the construction of appropriate wavelet bases on surfaces [10,13,14,23,28,31].In order to achieve the optimal complexity of the wavelet Galerkin scheme, wavelet bases are required with a sufficiently large number of vanishing moments.Our realization is based on biorthogonal spline wavelets derived from the multiresolution developed in [6].These wavelets are advantageous since the regularity of the duals is known [34].Moreover, the duals are compactly supported which preserves the linear complexity of the fast wavelet transform also for its inverse.This is an important task for the coupling of FEM and BEM, cf.[20,24,25].Additionally, in view of the discretization of operators of positive order, for instance, the hypersingular operator, globally continuous wavelets are available [2,7,13,23].
It turns out that the efficient computation of the relevant matrix coefficients is an important task for the successful application of the wavelet Galerkin method [23,29,31].We present a fully discrete Galerkin scheme based on numerical quadrature.Supposing that the given manifold is piecewise analytic we can use a hp-quadrature scheme [23,31,32] in combination with exponentially convergent quadrature rules.This yields an algorithm with asymptotically linear complexity without compromising the accuracy of the Galerkin scheme.
In view of nonsmooth geometries or singularities of the solution, a modern method should work adaptively.The most adaptive methods realize the adaptivity in the energy norm such that the super-convergence of the underlying Galerkin method is generally not realized.Wavelet bases offer the possibility to measure a wide range of Sobolev norms.In particular, adaptive schemes working with the optimal negative norm can be performed.Based on our actual approach we present these developments and provide numerical results which demonstrate the potential of our algorithm.
The outline of the present paper is as follows.First, in Section 2, we introduce the class of problems under consideration.Then, in Section 3 we provide suitable wavelet bases on manifolds.With such bases at hand we are able to introduce the fully discrete wavelet Galerkin scheme in Section 4. We survey on several practical issues like setting up the compression pattern, assembling the system matrix and preconditioning.In Section 5 we present an adaptive extension of the wavelet Galerkin scheme.In Section 6 we present numerical results in order to demonstrate our algorithm.

Problem Formulation and Preliminaries
For the numerical treatment of a boundary integral equation we need a discretization method which ends up with a sufficiently accurate finite-dimensional approximation of the given operator.At first we consider a general setting for the boundary element method.Next, a short description of the representation of the geometry on a computer is given.Then, we formulate the requirements to the kernel functions under consideration.

Boundary Integral Equations
We consider a boundary integral equation on the closed boundary surface Γ of a (n + 1)dimensional domain Ω where the boundary integral operator is assumed to be an operator of order 2q, that is The kernel functions under consideration are supposed to be smooth as functions in the variables x, y, apart from the diagonal {(x, y) ∈ Γ × Γ : x = y} and may have a singularity on the diagonal.Such kernel functions arise, for instance, by applying a boundary integral formulation to a second order elliptic problem.In general, they decay like a negative power of the distance of the arguments which depends on the spatial dimension n and the order 2q of the operator.
For the present purpose, we assume that the boundary Γ ∈ R n+1 is given as a parametric surface consisting of smooth patches, see Subsection 2.2 for details.The number of different mappings, which is the number of surface patches, will be denoted by M .The surface representation is in contrast to the usual approximation of the surface by panels.It has the advantage that the rate of convergence is not limited by this approximation.Moreover, such a representation is required for the adaptive solution of (1).Notice that technical surfaces generated by CAD tools are represented in this form.Of course, this fact makes the use of numerical integration indispensable for the computation of the system matrices.
The properties of the class of kernel functions k(x, y) which are under consideration will be outlined in Subsection 2.3.

Parametric Representation of Manifolds
Let denote the unit n-cube, i.e. = [0, 1] n .The manifold Γ ∈ R n+1 is partitioned into a finite number of patches where each γ i : → Γ i defines a diffeomorphism of onto Γ i .The intersection Γ i ∩ Γ i , i = i , of the patches Γ i and Γ i is supposed to be either ∅ or a lower dimensional face.
A mesh of level j on Γ is induced by dyadic subdivisions of depth j of the unit cube into 2 nj cubes C j,k ⊆ , where k = (k 1 , . . ., k n ) with 0 ≤ k m < 2 j .This generates 2 nj M elements Γ i,j,k := γ i (C j,k ) ⊆ Γ i , i = 1, . . ., M .
In order to ensure that the collection of elements {Γ i,j,k } on the level j forms a regular mesh on Γ, the parametric representation is subjected to the following matching condition: For each x ∈ Γ i ∩ Γ i exists a bijective, affine mapping Ξ : The first fundamental tensor of differential geometry is defined by the matrix Since γ i is supposed to be a diffeomorphism, the matrix K i (s) is symmetric and positive definite.The canonical inner product in L 2 (Γ) is given by The corresponding Sobolev spaces are denoted by H s (Γ), endowed with the norms • H s (Γ) , where for s < 0 it is understood that H s (Γ) = H −s (Γ) .Of course, depending on the global smoothness of the surface, the range of permitted s ∈ R is limited to s ∈ (−s Γ , s Γ ).

Kernel Functions and their Properties
We can now specify the kernel functions.To this end, we denote by α = (α 1 , . . ., α n ) and β = (β 1 , . . ., β n ) multi-indices of dimension n and define |α| := α 1 + . . .+ α n .Moreover, we denote by k i,i (s, t) the transported kernel functions, that is Definition 2.1.A kernel k(x, y) is called standard kernel of the order 2q, if the partial derivatives of the transported kernel functions k i,i (s, t), 1 ≤ i, i ≤ M , are bounded by We emphasize that this definition requires patchwise smoothness but not global smoothness of the geometry.The surface itself needs to be only Lipschitz.Generally, under this assumption, the kernel of a boundary integral operator A of order 2q is standard of order 2q.Hence, we may assume this property in the sequel.

Wavelets and Multiresolution Analysis
The nested trial spaces V j ⊂ V j+1 that we shall employ in the Galerkin scheme are spanned by so called single-scale bases Φ j = {φ j,k : k ∈ ∆ j } whose elements are normalized in L 2 (Γ) and whose compact supports scale like diam supp φ j,k ∼ 2 −j .
Associated with these collections and the choice of the wavelets are dual bases Φ j = { φ j,k : k ∈ ∆ j }, i.e., one has φ j,k , φ j,k = δ k,k , k, k ∈ ∆ j .For the current type of boundary surfaces Γ the Φ j , Φ j are generated by constructing first dual pairs of single-scale bases on the interval [0, 1], using B-splines for the primal bases and the dual components from [6] adapted to the interval [12].Tensor products yield corresponding dual pairs on .Using the parametric liftings γ i and gluing across patch boundaries leads to globally continuous single-scale bases Φ j , Φ j on Γ, [2,7,14,23].For B-splines of order d and duals of order for all s ≤ t ≤ d, s < γ, and likewise for the duals.It is known that the respective regularity indices γ, γ (inside each patch) satisfy γ = d − 1/2 while γ > 0 is known to increase proportionally to d.
Given the single-scale bases Φ j , Φ j , one can construct now biorthogonal complement bases [2,7,13,14] and [23] for particularly useful local representations of important construction ingredients.We suppose these complement bases normalized in L 2 (Γ).
A biorthogonal or dual pair of wavelet bases is now obtained by taking the coarse single-scale basis and the union of the complement bases where we have set for convenience Ψ j 0 −1 := Φ j 0 , Ψ j 0 −1 := Φ j 0 .Of course, in the infinite dimensional case the notion of basis has to be made more specific.The key feature of the wavelet basis is now the fact that Ψ, Ψ are actually Riesz bases in L 2 (Γ).In particular, the following norm estimate holds From biorthogonality and the fact that the dual spaces V j have the approximation order d one infers vanishing moments or the cancellation property of the primal wavelets Here The fact, that the concept of biorthogonality allows us to choose the order d of vanishing moments higher than the approximation order d, is essential for deriving optimal compression strategies that could not be realized by orthonormal bases.Appropriate piecewise constant and bilinear wavelets are depicted in Fig. 2.

The Wavelet Galerkin Scheme
This section presents a fully discrete wavelet Galerkin scheme for boundary integral equations.In the first subsection we discretize the given boundary integral equation.In Subsection 4.2 we introduce the a-priori matrix compression which reduces the relevant matrix coefficients to an asymptotically linear number.Then, in Subsection 4.3 and Subsection 4.4 we point out the computation of the compressed matrix.Next, in Subsection 4.5 we introduce an a-posteriori compression which reduces again the number of matrix coefficients.The last subsection is dedicated to the preconditioning of system matrices which arise from boundary integral operators of nonzero order.
In the sequel, the collection Ψ J with a capital J denotes the finite wavelet basis in the space V J , i.e., Ψ J := J−1 j=j 0 −1 Ψ j .Further, N J := dim V J ∼ 2 Jn indicates the number of unknowns.

Discretization
The variational formulation of the given boundary integral equation ( 1) reads For the Galerkin scheme we replace the energy space H q (Γ) in the variational formulation ( 6) by the finite dimensional spaces V J introduced in the previous section.Then, we arrive at the problem Equivalently, due to the finite dimension of V J , the ansatz ρ J = Ψ J ρ ψ J together with , yields the wavelet Galerkin scheme The system matrix A ψ J is quasi-sparse and might be compressed to O(N J ) nonzero matrix coefficients if the wavelets have a sufficiently large number of vanishing moments.The remainder of this paper is devoted to the efficient computation of the system matrix.
Remark 4.1.Replacing the wavelet basis Ψ J by the single-scale basis Φ J yields the traditional single-scale Galerkin scheme and ρ J = Φ J ρ φ J .This scheme is related to the wavelet Galerkin scheme by , where T J denotes the wavelet transform.The system matrix A φ J is densely populated.Therefore, the costs of solving a given boundary integral equation traditionally in the singlescale basis is at least O(N 2 J ).

A-priori compression
In a first compression step all matrix entries, for which the distance of the supports of the corresponding trial and test functions are larger than a level depending cut-off parameter B j,j , are set to zero.In the second compression step also some of those matrix entries are neglected, for which the corresponding trial and test functions have overlapping supports.First, we introduce the abbreviation Note that Ω j,k denotes the convex hull to the support of ψ j,k while Ω j,k denotes the socalled singular support of ψ j,k , i.e., those points where ψ j,k is not smooth.In accordance with [9,31] we have the following theorem.
Theorem 4.2.Let Ω j,k and Ω j,k be given as in (8) and define the compressed system matrix A ψ J , corresponding to the boundary integral operator A, by Fixing a > 1 and d < δ < d + 2q, the cut-off parameters B j,j and B j,j are set as .
Then, the error estimate Hence, in accordance with [35], the optimal order of convergence of the underlying Galerkin scheme is not compromised.The resulting structure of the compressed matrix is figuratively called finger structure, cf.Fig. 3.
The next theorem shows that the over-all complexity of assembling the compressed system matrix is O(N J ) even if each coefficient is weighted by a logarithmical penalty term [23].We mention that the choice α = 0 proves that the a-priori compression yields O(N J ) relevant matrix coefficients in the compressed system matrix.Theorem 4.3.Let the system matrix A ψ J be compressed according to Theorem 4.2.The complexity of computing this compressed matrix is O(N J ) provided that for some α ≥ 0 at most O J − j+j 2 α operations are spent on the approximate calculation of the nonvanishing coefficients (Aψ j ,k , ψ j,k ) L 2 (Γ) .

Setting up the compression pattern
Checking the distance criterion from Theorem 4.2 for each matrix coefficient, in order to assemble the compressed matrix, would require O(N 2 J ) function calls.To realize linear complexity, we exploit the underlying tree structure with respect to the supports of the wavelets, to predict negligible matrix coefficients.We will call a wavelet ψ j+1,son a son of ψ j,father if Ω j+1,son ⊆ Ω j,father .The following observation is an immediate consequence of the relations B j,j ≥ B j+1,j ≥ B j+1,j+1 , and B j,j ≥ B j+1,j for j > j .
With the aid of this lemma we have to check the distance criteria only for coefficients which stem from subdivisions of calculated coefficients on a coarser level.Therefore, the resulting procedure of checking the distance criteria is still of linear complexity.

Computation of Matrix Coefficients
Of course, the significant matrix coefficients retained by the compression strategy can generally neither be determined analytically nor be computed exactly.Therefore we have to approximate the matrix coefficients by quadrature rules.This causes an additional error which has to be controlled with regard to our overall objective of realizing asymptotically optimal accuracy while preserving efficiency.Theorem 4.3 describes the maximal allowed computational expenses for the computation of the individual matrix coefficients so as to realize still overall linear complexity.It is an immediate consequence of the fact that we require only a level dependent precision of quadrature, cf.[23,31].Theorem 4.5.Let the error of quadrature for computing the relevant matrix coefficients (Aψ j ,k , ψ j,k ) L 2 (Γ) be bounded by the level dependent accuracy with some d > d and δ ∈ (d, d + r) from Theorem 4.2.Then, the Galerkin scheme is stable and converges with the optimal order (9).
From (10) we conclude that the entries on the coarse grids have to be computed with the full accuracy while the entries on the finer grids are allowed to have less accuracy.Unfortunately, the domains of integration are very large on coarser scales.Since we employ primal multiresolution spaces V j based on piecewise polynomials, the numerical integration can be reduced to the computation of the interaction of polynomial shape functions on certain elements.Consequently, we have only to deal with integrals of the form with p l denoting the polynomial shape functions and the transported kernel function (3).This is quite similar to the traditional Galerkin discretization.The main difference is that in the wavelet approach the elements may appear on different levels due to the multilevel hierarchy of wavelet bases.Difficulties arise if the domains of integration are very close together relatively to their size.We have to apply numerical integration with some care in order to keep the number of evaluations of the kernel function at the quadrature nodes moderate and to fulfill the requirements of Theorem 4.3.The necessary accuracy can be achieved within the allowed expenses if we employ an exponentially convergent quadrature method.
In [23,31] a geometrically graded subdivision of meshes is proposed in combination with varying polynomial degrees of approximation in the integration rules, cf.Fig. 4. Exponential convergence is shown for boundary integral operators under the assumption that the underlying manifolds are piecewise analytic.It is shown in [23] that the combination of tensor product Gauß-Legendre quadrature rules with such a hp-quadrature scheme leads to a number of quadrature points satisfying the assumptions of Theorem 4.3 with α = 2n.Since the proofs are rather technical we refer to [32,29,31,23] for further details.For that result to be valid we need a slightly stronger assumption on our manifold Γ which should be piecewise analytic.Moreover, the kernels of the operators should satisfy the following condition.
Assumption 4.6.The kernel function k(x, y) is analytically standard of order 2q, that is, the partial derivatives of the transported kernel function (3) are uniformly bounded by with some q > 0.
Since the kernel function has a singularity on the diagonal we are still confronted with singular integrals if the domains of integration live on the same level and have some points in common.This happens if the underlying elements are identical or share a common edge or vertex.When we do not deal with weakly singular integral operators, the operators can be regularized, e.g. by partial integration [27].So we end up with weakly singular integrals.Such weakly singular integrals can be treated by the so-called Duffy-trick [16,30].In this way the singular integrands are transformed into analytical ones.

A-posteriori compression
If the coefficients of the compressed system matrix A ψ J have been computed, we may apply an a-posteriori compression by setting all coefficients to zero, which are smaller than a level depending threshold, cf.[9,23].That way, a matrix A ψ J is obtained which has less nonzero coefficients than the matrix A ψ J .Clearly, this does not affect the calculation of the matrix coefficients but accelerates the time for the iterative solution.Moreover, the requirement to the memory is reduced if the system matrix has to be stored.For instance, this is advantageous for the coupling of FEM and BEM, cf.[20,24,25].To our experiences this procedure reduces the number of nonzero coefficients by a factor 2-5.
Theorem 4.7.We define the a-posteriori compression by Herein, the level dependent threshold ε j,j is chosen as in (10) with some d > d and δ ∈ (d, d + r) from Theorem 4.2.Then, the optimal order of convergence (9) of the Galerkin scheme is not compromised.

Wavelet preconditioning
If A : H q (Γ) → H −q (Γ) is a boundary integral operator of nonzero order, the corresponding system matrix A ψ J is ill conditioned.In fact, there holds cond l 2 A ψ J ∼ 2 2J|q| .According to [8,11,31], the wavelet approach offers a simple diagonal preconditioner based on the norm equivalences.
Theorem 4.8.Let the diagonal matrix D r J defined by Then, if A : H q (Γ) → H −q (Γ) denotes a boundary integral operator of the order 2q with γ > −q, the diagonal matrix D 2q J defines a preconditioner to A ψ J , i.e., Remark 4.9.The coefficients on the main diagonal of A ψ J satisfy Therefore, the above preconditioning can be replaced by a diagonal scaling.In fact, the diagonal scaling improves and simplifies the wavelet preconditioning.
As the numerical results in [26] confirm, this preconditioning works well in the two dimensional case.However, in the three dimensions, the results are not satisfactory.One figures out of Fig. 5 the condition numbers of the stiffness matrices of the single layer operator on a square discretized by piecewise bilinears.We employed different constructions for wavelets with four vanishing moments (spanning identical spaces, see [23] for details).In spite of the preconditioning, the condition numbers with respect to the wavelets are not significantly better than with respect to the single-scale basis.We mention that the situation becomes even worse for operators defined on more complicated geometries.A slight modification of the wavelet preconditioner yields much better results [23].The simple trick is to combine the above preconditioner with the mass matrix which yields an operator based preconditioning.Theorem 4.10.We consider a boundary integral operator A : H q (Γ) → H −q (Γ) with corresponding Galerkin matrix A ψ J .Let D r J be defined as in (12) and denote the mass matrix.Then, if γ > −q, the matrix C 2q J = D q J B ψ J D q J defines a preconditioner to A ψ J , i.e.
This preconditioner decreases the condition numbers impressively, cf.Fig. 5. Let us remark that the condition depends on the underlying spaces but not on the chosen wavelet basis.To our experiences the condition reduces about the factor 10-100 compared to the preconditioner (12).

Adaptivity
Wavelet matrix compression can be viewed as a non-uniform approximation of the Schwartz kernel k(x, y) with respect to the typical singularity at x = y due to Defintion 2.1.If the domain has corners and edges, the solution itself admits singularities.In this case an adaptive scheme will reduce the number of unknowns drastically without compromising the overall accuracy.Adaptive methods for BEM have been considered by several authors, see e.g.[3,17,18,19,33] and the references therein.However we are not aware of any paper concerning the combination of adaptive BEM with recent fast methods for integral equations like e.g. the fast multipole method.Similar to the matrix compression, in the wavelet context a non-uniform or adaptive approximation is achieved by compression, i.e., deleting small wavelet coefficients.Therefore an adaptive approximation space V J , spanned by the wavelet bases with active indices, is clearly a subspace of the space V J .Consequently, we pursue to find a sequence of spaces where V j ⊆ V j , such that the Galerkin solution with respect to V j provides nearly the same accuracy as the Galerkin solution with respect to V j .
In this respect the combination of matrix compression with adaptivity via vector compression seems to be quite natural.Recent work of Cohen et al. [4,5] is directed to this strategy.These authors showed how to achieve a numerical solution ρ up to accuracy with N ∼ −1/s unknowns in O(N ) complexity.Their algorithm is based on an iterative method for the continuous equation expanded with respect to a wavelet basis.The application of the operator to a function is approximated by an appropriate finite matrixvector multiplication.From time to time the actual approximate solution is compressed by thresholding providing the best N -term approximation.Since the computation of the relevant matrix coefficients is by far the most expensive step in our algorithm, we cannot use this approach of [4,5] directly.At least we adopt the strategy of the best N -term approximation by wavelet compression combined with the notion of tree approximation as considered in [15].
To define appropriate wavelet basis on V j , we have to ensure that the support of a small wavelet does not intersect large elements.This appears if small and large elements are immediately adjacent.We call a mesh 1-graded, if the levels of adjacent elements differ at most by one.Likewise, a graded mesh is 2-graded, if the levels of all neighbours of each element differ at most by one, and so on, cf.Fig. 6.An exact definition of m-gradedness can be found in [15].The gradedness ensures that we find a tree structured (with respect to the supports of the associated wavelets) index set Λ j such that V j = span{ψ λ : λ ∈ Λ j } with |Λ j | = dim V j .Moreover, completing Λ j by the sons of all leaves, we obtain the index set Λ j, , which generates the trial space V j, = span{ψ λ : λ ∈ Λ j, } that arises from uniform subdivision of V j .We mention that the mesh has to be patchwise 1-graded in the case of the piecewise constant wavelets presented in Fig. 2. The piecewise bilinear wavelets require 2-gradedness, which has to be extended to global 2-gradedness if we consider them globally continuous.
The rate of convergence in the adaptive algorithm of [4,5] was considered in the energy norm.To our experience this rate is too low because it restricts in general the superconvergence of the Galerkin scheme.The highest order of convergence of the boundary element method is achieved with respect to the norm in H 2q−d (Γ), cf.(9).Since the number of vanishing moments is chosen such that d < d + 2q, we can estimate this norm by (4).If we consider the saturation assumption as for example in [18,19], the optimality of our method follows immediately.Assumption 5.1.Let V j denote an arbitrary m-graded trial space and let V j, be the trial space that arises from uniform subdivision of V j .For a given t ∈ [2q − d, γ) we assume that that there exists a constant θ < 1 such that the solutions ρ j with respect to V j and ρ j, with respect to V j, satisfy ρ − ρ j, t ≤ θ ρ − ρ j t .
Up to now we can compute the Galerkin solutions ρ j and ρ j, .Our problem reads now: find the smallest index set Λ j ⊆ Λ j+1 ⊆ Λ j, , such that the Galerkin solution ρ j+1 with respect to ψ Λ j+1 satisfies (14).At present we choose the canonical strategy and compute elementwise error portions by bunching the wavelets which correspond to the subdivision of an element of V j .This procedure is simple to implement and corresponds to that when using hierarchical error estimators, see [18,19].However, we use wavelet thresholding instead of a-posteriori error estimators.After sorting these error portions by their modulus, we increase the index set Λ j successively by activating the wavelets corresponding to the largest error portions until ( 14) is satisfied.Possibly the so constructed index set Λ j+1 has to be extended to ensure the m-gradedness of the new trial space V j+1 .
We are now in the position to formulate our adaptive algorithm, which is based on a nested iteration.initialization: V 0 := V 0 for j := 1 to J − 1 do begin compute the compressed system matrix for V j−1, compute the solutions ρ j−1 and ρ j−1, determine V j such that (14) holds end compute the compressed system matrix for V J−1, compute the final solution ρ J := ρ J−1,
We consider the gearwheel presented in Fig. 7 as domain Ω.Choosing the harmonical function u(x) = [4, 2, 1]x x −3 and setting f = u| Γ , the Dirichlet problem has the unique solution u.
We make use of the indirect formulation involving the single layer operator V : which gives a Fredholm integral equation of the first kind Vρ = f on Γ.
We solve this boundary integral equation by the traditional Galerkin scheme using piecewise constant ansatz functions and by the nonadaptive and adaptive wavelet Galerkin scheme Figure 7: The mesh of a gearwheel parametrized by 172 patches after two subdivision steps.
using piecewise constant wavelets with three vanishing moments.The setting for the adaptive scheme is t = −2 and = 1/3 First, in Tab. 1 we tabulate the maximum norm of the absolute errors of u J and the overall cpu-times T cpu (measured in seconds).The optimal order of convergence is cubic, i.e. ε ∼ h 3 J ∼ dim V J −3/2 , but we expect approximately only a quadratic order ε ∼ h 2 J due to reentrant edges and vertices of the geometry.Notice that level 3 and 5 is no more computable with the traditional Galerkin scheme and the (nonadaptive) wavelet Galerkin scheme, respectively.All schemes provide a nearly identical accuracy as far as comparable, but the adaptive one does it with essentially less unknowns and, hence, less memory requirement and cpu-time.The improvement of the adaptive over the nonadaptive wavelet Galerkin scheme is similar to the improvement of the (nonadaptive) wavelet Galerkin scheme over the traditional Galerkin scheme.Fig. 8 visualizes the compression rates of the (nonadaptive) wavelet Galerkin scheme.We plot the number of nonzero coefficients in percent.For 74240 unknowns the matrix compression yields only 1.8 % and 0.49 % relevant matrix coefficients after the a-priori and a-posteriori compression, respectively.
Next, in Fig. 9 we plot the ratio dim V J / dim V J corresponding to the adaptive wavelet Galerkin scheme.As the dashed curve indicates, we conclude an approximate increasement dim V J+1 = 2.5 dim V J of the adaptive spaces instead of dim V J+1 = 4 dim V J of the full wavelet spaces.The adaptive scheme should converge with the order ε ∼ h 2 J , where h J is the mesh size of the finest grid.Relating the accuracy to the number of unknowns N J = dim V J , we observe ε ∼ N −3/2 J in accordance with Tab. 1.This is exactly the optimal rate of convergence because smooth solutions can be approximated only by this rate ε ∼ N −3/2 J .Such a behaviour is predicted by the papers [4,5], however the present rate is better than the rate with respect to the energy norm as in [4,5].

Figure 1 :
Figure 1: The parametrization of the unit sphere is obtained by projecting it onto the cube [−1, 1] 3 , which yields six patches (left).On the right hand side one figures out the partition on the level j = 4.

wavelet of type i wavelet of type ii wavelet of type iiiFigure 2 :
Figure 2: (Interior) piecewise constant/bilinear wavelets with three/four vanishing moments.

Figure 3 :
Figure 3: The finger structure of the compressed system matrix computed with respect to the two dimensional (left) and the three dimensional (right) unit spheres.

Figure 4 :
Figure 4: Adaptive subdivision of the domains of integration.
single−scale basis diagonal scaling: tensor product wavelets diagonal scaling: simplified tensor product wavelets diagonal scaling: wavelets optimized w.r.t. the supports modified preconditioner

Figure 5 :
Figure 5: The l 2 -condition numbers with respect to the single layer operator on the unit square and piecewise bilinear wavelets with four vanishing moments.

Figure 8 :
Figure 8: Compression rates of the system matrix of wavelet Galerkin schemes.
unknowns in the adaptive scheme Level J number of active unknowns in per cent

Figure 9 :
Figure 9: Number of unknowns of the adaptive wavelet Galerkin scheme related to the nonadaptive one.

Table 1 :
. The approximate potentials u J = [(Vρ J )(x i )] are calculated in many points x i distributed inside the gearwheel.The exact potential is denoted by u = [u(x i )].The computations are performed by a standard personal computer with 1 Gigabyte main memory.Accuracy and cpu-times.