• No results found

Using wavelets in the dual reciprocity method

N/A
N/A
Protected

Academic year: 2021

Share "Using wavelets in the dual reciprocity method"

Copied!
23
0
0

Bezig met laden.... (Bekijk nu de volledige tekst)

Hele tekst

(1)

Using wavelets in the dual reciprocity method

Citation for published version (APA):

Morsche, ter, H. G., Mattheij, R. M. M., & Wang, K. (1999). Using wavelets in the dual reciprocity method. (RANA : reports on applied and numerical analysis; Vol. 9924). Technische Universiteit Eindhoven.

Document status and date: Published: 01/01/1999 Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers) Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

(2)

Using Wavelets in the Dual Reciprocity Method

H.G. ter Morsche, R.M.M. Mattheij & K.Wang

EMail: morscheh@win.tue.nl

Abstract

Radial basis functions, in especial thin plate splines or conic splines seem to be very useful in the Dual Reciprocity Method (DRM) for solving partial dierential equations by means of BEM. In DRM the radial basis function have been used for computing a particular solution for an inhomogenuous partial dierential. Two main steps are important in applying DRM. First, a functionf must be approximated by shifted radial

basis functions. Second, a particular solution must be found in an easy way for the radial basis function itself. Despite the fact that radial basis functions may have good approximation properties, the approximation scheme's involved are far from local. The consequence is that many coecients are needed to represent functions by series of shifted radial basis functions, which causes a lot of computational eorts. Wavelets have the potential to overcome this problem. However, for wavelets the problem of nding a particular solution looks more complicated. In this chapter, it will be shown that for the so-called hexagonal wavelets in two-dimensions and the Poisson equation, it is still possible to derive a particular solution by doing elementary function evaluations.

1 Introduction

Consider the Poisson equation in two dimensions:

4u(x) = f(x) (x2): Here x= (x1x2) 2IR 2,  = @ 2 @2 x 1 + @2 @2 x 2

is the 2D Laplacian, f is a presribed bivariate function and IR

2 is a bounded simply connected domain in IR2. To

nd the functionu(xy), which satis es the Poisson equation on  subject to the Dirichlet condition

u(x) =g(x) (x2;):

where ; is the closed piecewise smooth boundary of  and g a prescribed continuous function, by means of BEM methods, one rst convert this problem into an integral as follows:

LetG(x) be the fundamental function corresponding to the Laplacian4. So, 4G(x) =(x), where (x) is the bivariate Dirac delta function. One has

(3)

G(x) = 12 logjxj (1) wherejxj= q x2 1+x 2

2. It follows from Green's second formula that for points

x in the interior of  the following holds.

u(x); Z  Gx(y)f(y)d = Z ; u(y)@G@n dx(y) ;; Z ; Gx(y)@u@n d(y) ;: (2)

whereGx(y) stands forG(x;y) and wheren is the outer normal of .

Whenx belongs to the boundary ;, then depending on the local smoothness of ; at x, the value u(x) in the foregoing equation must be replaced byc(x)u(x) for some factor c(x). For instance if ; isC1 at x, then c(x) = 1=2. So if x is a

boundary point, we have

c(x)u(x); Z  Gx(y)f(y)d = Z ; u(y)@G@n dx(y) ;; Z ; Gx(y)@u@n d(y) ;

which in case of our Dirichlet condition yields an integral equation for the normal derivative of u on the boundary. If this equation has been solved, then we may apply formula (2) to obtain the solutionuin the interior of . Evidently, the same strategy can be used for the Neumann problem (the Dirichlet boundary condition

u(x) = g(x) is replaced by the Neumann boundary condition @u(x)

@n = g(x))

and even for mixed versions of the Dirichlet and Neumann conditions involving dierent parts of the boundary (cf. 5]).

The boundary integral equation (3) contains a domain integral, this is the integral over . If a function up can be found such that 4up = f, then before

applying numerical methods to solve the integral equation, we rst reduces the domain integral to boundary integrals by again using Green's second formula. This is, what P.W. Partridge e.a. call in their book (5]) the Dual Reciprocity Method (DRM). The DRM avoids computing domain integrals. For our Dirichlet problem this leads to:

Z  Gx(y)f(y)d = Z  Gx(y)4up(y)d = c(x)up(x) + Z ; Gx(y)@u@n dp(y) ;; Z ; up(y)@G@n dx(y) ;:

So, the integral equation that must be solved has the form:

c(x)g(x) = c(x)up(x) + Z ; Gx(y)@u@n dp(y) ;; Z ; up(y)@G@n dx(y) ; +Z ; g(y)@G@n dx(y) ;; Z ; Gx(y)@u@n d(y) ;: (3)

(4)

To solve this boundary integral equation by means of nite boundary ele-ments, one normally needs the values of up(x) and the gradient rup(x) at a

certain set of boundary points. If elementary function evaluations for up(x) and

rup(x) are available, then these values can be computed in a direct way.

In practice the functionf will be approximated by a superpostionf =P

kfk

of functions fk such that for each functionfk an analytic expression for a

corre-sponding particular solution can be found straightforwardly. Popular choices are radial basis functions fk(x) = (jx;xkj) , where (r) is a function de ned on

IR+ and (x

k) a collection of points, which can be chosen in some way. Since the

Laplacian is translation invariant, we only have to look for a particular solution of the equation u =(jxj), which may also be considered as a function of jxj

only. Finding a particular solution is then solving an ordinary 1D dierential equation. In 1995, M.A. Golberg (4]) gave a survey of the post-1990 literature on the numerical evaluation of particular solutions in the BEM until that time. With respect to the Poisson equation radial basis functions, in especial thin plate splines or conic splines (cf. 6]), seem to be very useful. The disadvantage is the complexity of the associated approximation problem, since the radial basis functions are in general not compactly supported.

In this chapter we replace the radial basis functions by compactly supported 2D-wavelets, called hexagonal wavelets, which are radial symmetric in some sense. Then we may pro t from the fact that in a wavelet expansion of a function

f, a relatively small number of the wavelet coecients are needed to represent the function f. For our applications this means, that by clipping the wavelet coecients of f with are small compared to a certain treshold value, a vast amount of coecients in the wavelet expansion of f can be set to zero, which reduces the computational eort for setting up the boundary integral equation.

In this chapter, we are concerned with the problem of nding a particular solution of the inhomogenous Poisson equation when these hexagonal wavelets are used. It turns out that a particular solution can be expressed in terms of third order dierences of a function at the vertices of a cube in IR3. Section 2

contains an overwiev of the theory of wavelets in two dimensions, which is needed to understand the hexagonal wavelets presented in Section 3. We also give in Section 3 some examples of the consequences of tresholding wavelet coecients in the wavelet expansion of some functions. Finally at Section 4, particular solutions will be derived corresponding to the hexagonal wavelets and a numerical test will be carried out.

2 Wavelets preliminaries

One of the applications of wavelet decompositions in numerical analysis is the sparsi cation of full matrices which may occur in the discretisation of integral

(5)

equations. It is based on the property that a relative high percentage (depending on the choice of the wavelets) of the wavelet coecients in the wavelet expansion of suciently smooth functionsf 2L

2(IR2) is closed to zero. In order to

under-stand this property of wavelets, which is important for our contribution to the dual reciprocity method, this section will contain an introductory overview of the theory of 2D-wavelets. For (a lot) more information the reader is referred to text books on wavelets. We mention a few of them here: Daubechies (3]), Chui (1]), Sidney Burrus e.a (7]) and Strang (8]).

For bivariate functions a wavelet expansion is normally based on three so-called mother wavelets 1 2 and 3, from which by translations and dyadic

scaling the wavelets are derived. The corresponding wavelet expansion usually has the form:

f(x) = 1 X j=;1 3 X l=1 X n2ZZ 2 dlnj lnj(x) (4) where the wavelet coecients dlnj have nitel2-norm, i.e.,

3 X l=1 1 X j=;1 X n2ZZ 2 jdlnjj 2 < 1:

In this expansion the functions lnj, called wavelets, are de ned as follows:

lnj(x) = 2j l(2jx;n 1e1 ;n 2e2) (x 2IR 2): (5) Here n= (n1n2) 2 ZZ 2 and e

1 and e2 are two independent vectors in IR 2,

which may dier from the unit vectors (10) and (01). Note thatk lnjk=k lk,

where kkdenotes the L

2-norm in the function spaceL2(IR2).

Apparently, the function f is written as a superposition of translated and scaled versions of the three so-called mother wavelets. It is required that

ZZ

IR 2

l(x)dx= 0 (l= 123): (6)

So, the mean value of all the wavelets lnj are equal to zero.

The wavelet expansion of f can also be seen as a decomposition of f into functions fj corresponding to dierent space scales. To be more precise, we may

write

f = 1 X

j=;1

(6)

where fj(x) = P 3

i=1 P

n2ZZ

2dlnj lnj(x) \represents" f in the scale space

Wj IR

2 spanned by the wavelets f lnj  l = 123 n 2 ZZ 2 g, where j 2 ZZ

is xed. The spaces Wj are closed linear subspaces of L2(IR

2) and, in fact,

di-lated versions of the space W0:

f(x)2W

0 if and only if f(2

jx)2Wj:

Moreover one has: (i) WjT Wj =f0g (j6=i), (ii) L2(IR 2) = +W ;2+W;1+W0+W1+W2+ .

Here the +-sign refers to direct summation of spaces. In case the basis

f lnj  l= 123 j 2Z n2ZZ 2 g is an orthogonal basis inL 2(IR 2), the wavelets

are called orthogonal wavelets. The well-known Daubechies family of wavelets delivers examples of orthogonal wavelets having nite supports. However, or-thogonal wavelets having compact supports lack symmetry properties. A way to overcome this problem is to leave the orthonormality property and to consider the more general class of the so-called bi-orthogonal wavelets. In the context of bi-orthogonal wavelets, one is faced with two wavelet bases of L2(IR2), the

-basis and the ~ basis, which are mutually related in the following way: ( lnj ~kmi) =lknmji:

For orthonormal wavelet bases, these two bases coincide. It follows straight forward from the de nition of bi-orthogonal wavelet bases that the wavelet coef- cients dlnj in (4) are given by the innerproducts:

dlnj= (f ~lnj):

For stability purposes it is important that the wavelet bases is a so-called Riesz-basis. This means that positive constants A andB exist such that for the wavelet coecients in (4) the following inequalities hold:

A 1 X j=;1 3 X l=1 X n2ZZ 2 jdlnjj 2 kfk 2 B 1 X j=;1 3 X l=1 X n2ZZ 2 jdlnjj 2

If the -basis is a Riesz basis then also the ~ basis is a Riesz- basis. The constants A and B are then replaced by 1=B and 1=A respectively. We now return to our scale spaces Wj.

Note that Wj is a direct sum Wj =W1

j +W2 j +W3 j of the spaces W1 j, W2 j and W3

j which are dilated versions of W1 0 W

2

0 and W 3

0 respectively. The space

Wl

0 is spanned by the wavelets

f l(x;n 1e1 ;n 2e2)  n 2ZZ 2 g. If we add (by a

(7)

direct sum) all the scaled spaces Wj up to the scalej =k;1, then one obtains

a sequence of spaces (Vk) (k 2ZZ), given by

Vk= k

;1 X

j=;1

Wj:

These spaces satisfy the properties: (i) V ;1 V 0 V 1 ,

(ii) f(x)2Vk if and only if f(2x)2Vk +1, (iii) \ k2ZZ Vk =fog, (iv) k2ZZ Vk =L2(IR2).

For the most common wavelets the following property holds: (v) A functionexists such that

f(x;n 1e1 ;n 2e2)  n= (n 1n2) 2ZZ 2 g is a Riesz-basis ofV0.

The sequence of spaces (Vk) is called a multiresolution analysis (MRA) of

L2(IR

2) and its scaling function. Often, the scaling function is normalized by

the condition

ZZ IR

2

(x)dx= 1: (7)

However, we prefer to have scaling functions, like the B-splines in the ap-proximation theory, being a partition of unity. This means:

X n1n2 (x;n 1e1 ;n 2e2) = 1 (x 2IR 2): (8)

The wavelets, we will apply in this chapter have a nonnegative scaling func-tion, so one has

ZZ IR

2

(x)dx >0: (9)

In general, the construction of wavelets initiates by setting up a multi-resolution analysis to which an appropriate scaling function is associated. The dual wavelets stem from the dual MRA, which is described by a sequence of spaces ~Vk and associated scaling function ~, which has the property:

(n0~m0) =nm:

(8)

-

Wi ?W~j(i6=j)

-

Vk =Vk;1+W 1 k;1+W 2 k;1+W 3 k;1

-

V~k = ~Vk;1+ ~W 1 k;1+ ~W 2 k;1+W 3 k;1

-

Wlk is spanned by f lnk  n 2ZZ 2 g,

-

W~lk is spanned by f ~lnk  n 2ZZ 2 g,

-

Vk is spanned byfnk  n 2ZZ 2 g,

-

V~lk is spanned by f~nk  n 2ZZ 2 g,

>From this list a lot of mutual relations between the wavelets and the scaling functions can be obtained. These relations will show how close the wavelets are related to the theory of lter banks.

Since 2V 0

V

1, there exists a sequence of numbers (q 0 n) such that (x) = X n2ZZ 2 q0 nn1(x): (10)

It can be shown that for scaling functions having nite supports, only a nite number of elements in the sequence (q0

n) diers from zero. In such a situation a

sequence is said to be nite. We like to have nitely supported scaling functions and wavelets, so we have to assume that the sequence (q0

n) is nite.

Applying the 2D Fourier-transform ( f(x) ! f^(!) = R R IR

2

f(x)e;i!xdx) to

the foregoing relation, one has: ^ (!) = 1 2 X n2ZZ 2 q0 ne;i(n 1! e 1 +n 2! e 2 )=2^(! 2): (11)

Here the dotstands for the inner-product in IR

2. Similarly, since lhas nite

support and l 2 Wl 0

 V

1, there exist three nite sequences (qln) (l = 123)

such that

l(x) = X

n2ZZ 2

qlnn1(x): (12)

Again using the 2D Fourier-transform this yields ^ l(!) = 1 2 X n2ZZ 2 qlne;i(n 1! e 1 +n 2! e 2 )=2^l(! 2):

(9)

We introduce now the functions Ml(!) as follows: Ml(!) = 12 X n2ZZ 2 qlne;i(n 1! e 1 +n 2! e 2 ) (l= 0123): (13)

The functionsMl(!) are 2-periodic with respect to the dual pair of vectors

~ e1 and ~e2 de ned by e1 e 1 = 1e~1 e 2= 0 ~ e2e1 = 0e~2 e 2= 1: (14) So Ml(!+ 2e~1) =Ml(!+ 2e~2) =Ml(!)(! 2IR

2). Moreover, one has

^ (!) = M0( ! 2 )^(!2 ) (15) ^ l(!) = Ml(! 2 )^(!2 ) (l= 123): (16) Similar relations hold for the dual wavelets and scaling function. To be complete we list them here:

~ (x) = X n2ZZ 2 ~ q0 n~1n(x) (17) ~ l(x) = X n2ZZ 2 ~ qln~1n(x) (l= 123) (18) ^~ (!) = ~M0( ! 2 )^~(!2 ) (19) ^~ l(!) = ~Ml(!2 )^~(!2 ) (l= 123): (20)

Also the sequences (~qln) (l= 0123) are assumed to be nite. It follows form (6) and (9) that

^

l(0) = 0 (l= 123)

^

(0)6= 0:

It is also required that ^~

l(0) = 0 (l= 123)

^~

(10)

So, one has

Ml(0) = ~Ml(0) = 0 (l= 123)

M0(0) = ~M0(0) = 1:

The duality relations may also be expressed by means of Fourier-transforms. By setting 0= (00)1 =e~12 =~e2 (cf. (14)) and3 =(~e1+ ~e2), it can

be shown that 3 X l=0 Mj(!+l) ~Mk(!+l) =jk (!2IR 2j 6 =k): (21) The functions Ml(!) and ~Ml(!) play an important role in the description

of wavelets with help of lter banks. The lter banks are used to decompose a function in wavelets and also to reconstruct a function from its decomposition.

Decomposition and reconstruction

In practice, the wavelet coecients are seldom computed by taking the inner products with the dual wavelets. If a functionf is given, then rst f will be ap-proximated by a functionf0

2V

0 say, by means of an approximation scheme for

instance by interpolation. The space V0 must have the property that it contains

a \good" set of candidates to approximate the given function f. In the next section, the space V0 will consist of functions which are continuous and

piece-wise linear on a regular mesh with small mesh size. By taking into account the approximation error, we replace f by f0. So, we have

f(x) f 0(x) = X n2ZZ 2 an0n0(x) = ;1 X j=;1 3 X l=1 X n2ZZ 2 dljn lnj(x) (x2IR 2): Since,V0 =V;1+W 1 ;1+W 2 ;1+W 3 ;1, we may write f0(x) = X n2ZZ 2 an;1n;1(x) + 3 X l=1 X n2ZZ 2 dln;1 ln;1(x): (22)

How the new approximation coecients (an;1) and the wavelet coecients

(dln;1) (also called detail coe cients in this context) at level

;1 depend on the

approximation coecients (a0

(11)

an;1 = X k2ZZ 2 ~ q0 k;2nak 0 dln;1 = X k2ZZ 2 ~ qlk;2nak 0 (l= 123):

Note the occurence of the coecients (~qln) (cf. 18) in these formulas. We also observe that the coecients at level ;1 are of the form

P

k2ZZ 2q~lk

;mak

0, where

m = 2n. So, rst one has to compute the convolution product of the sequences (an0) and (~q

l

;n), which can be considered as a linear lter proces (convolution)

with frequency response

X n2ZZ 2 ~ ql ;ne ;i(n 1! e1+n 2! e2)= 2M l(!):

Then the \output sequence ", (bn) say, is downsampled by restricting to the

elements (b2n). The lter proces + downsampling is schematically presented in

Figure 1 -2 ~M1(!) 2 ~M2(!) (an0) -(d1 n;1) (d2 n;1) ? ? -2 ~M0(!) 2 ~M3(!) -? ? ? -(d3 n;1) -(an;1) Figure 1: decomposition On the other hand, if the coecients (an;1), (dln

;1), are given then we may

(12)

c0 n0 = X k2ZZ 2 q0 n;2kak ;1 c1 n0 = X k2ZZ 2 q1 n;2kd 1 k;1 c2 n0 = X k2ZZ 2 q2 n;2kd 2 k;1 c3 n0 = X k2ZZ 2 q3 n;2kd 3 k;1 an0 = c 0 n:0+c 1 n0+c 2 n0+c 3 n0

Also, this can be considered as a linear lter proces, which is followed by summation. However, before the lter process is applied, an upsampling is needed by inserting zeros. Upsampling of a sequence ( n) gives a new sequence (n)

such that 2n= nfor alln 2ZZ

2. The reconstruction is schematically shown in

Figure 2. The frequency responses of the lters are given by 2M0(!), 2M1(!),

2M2(!), 2M3(!) respectively. -(d1 n;1) (d2 n;1) (d3 n;1) (an;1) 6 6 6 6 -2M3(!) 2M2(!) 2M1(!) 2M0(!) -+ -(an0) Figure 2: reconstruction

The decomposition lters decompose the function f0 in four functions the

new approximation function f;1 and at level

;1 the three detail functions g 1 ;1, g2 ;1 and g 3 ;1. So,f 0 =f;1+g 1 ;1+g 2 ;1+g 3

;1. We may continu, using the same

decomposition lters, to decompose the functionf;1 into the functionsf;2,g 1 ;2, g2 ;2, and g 3 ;2. Then we have: f 0=f;2+g 1 ;2+g 2 ;2+g 3 ;2+g 1 ;1+g 2 ;1+g 3 ;1. Of

(13)

course, we still may go on and decompose the function f;2 etc. Since, g

l

;1(l=

123) contains \information" of f0 on the nest available scale, a lot of detail

coecients in (dln;1) are expected to be close to zero. By choosing an appropriate

treshold value these coecients can be replaced by zero. This also may happen for the functions gl

;2, etc. In image processing, this proces is applied as an

eective method to compress images. In numerical analysis, this method can be used to reduce the number of parameters for the representation of a function, which generally will lead to less computations. An indication for choosing an appropriate can be the number of vanishing moments of the dual wavelets ~ l

(l= 123). This is the maximal numberm for which

Z IR

2

p(x) ~ l(x)dx= 0 (l= 123) (23)

for all polynomials p(x) of (total) degree at mostm. Since dln;1= Z IR 2 f(x) ln;1(x)dx

it follows from Taylors formula and the de nition of m that for functions f(x) wich are m+ 1 times continuous dierentiable the following inequality holds

jdln ;1 jC Z IR 2 jx;ajm +1 j ln ;1(x) jdx (24)

where ais any point in the support of ln;1.

The functions ~ l we shall apply, have nite support, so it follows from (24)

that jdln ;1

j is of the order of (diam(support( ln ;1))

m+3, where diam is the

di-ameter of the support of l. So if the support of ~ l is small, we may expect small

values for the numbersdln;1.

In the next section, we will focuss on speci c wavelets, called hexagonal wavelets. Moreover, an one level decomposition will be applied to the Gaussian function e;x 2 1 +x 2 2.

3 Hexagonal wavelets

The most frequently used bivariate wavelets have been constructed by means of tensorproduct of univariate wavelets. In this case, the starting point is a bi-orthogonal system of wavelets generated by a univariate motherwavelet (x1) (we

allow ourselves to use the same function symbols as in the two-dimensional case, only the number of arguments dier) its dual ~ (x1) and the univariate scaling

(14)

functions(x1) and ~(x1). Then, from these functions the bivariate wavelets and

scaling functions can be obtained as follows:

(x) =(x1)(x2)~(x) = ~(x1) ~(x2) 1(x) =(x 1) (x2) ~ 1(x) = ~(x 1) ~ (x2) 2(x) = (x 1)(x2) ~ 2(x) = ~ (x 1) ~(x2) 3(x) = (x 1) (x2) ~ 1(x) = ~ (x 1) ~ (x2): Here x= (x1x2).

In one dimension, one has two reconstruction ltersp

2M0(!1) and p

2M1(!1)

with two decomposition lters p

2 ~M0(!1) and p

2 ~M1(!1).

These lters generate the two dimensional lters using the formulas:

M0(!) =M0(!1)M0(!2)M~0(!) = ~M0(!1) ~M0(!2)

M1(!) =M0(!1)M1(!2)M~1(!) = ~M0(!1) ~M1(!2)

M2(!) =M1(!1)M0(!2)M~2(!) = ~M1(!1) ~M0(!2)

M3(!) =M1(!1)M0(!2)M~3(!) = ~M1(!1) ~M2(!2)

Here != (!1!2).

Apparently, in the tensor case a two dimesional lter is a cascade of two one dimensional lters. In one dimension, the ltersM0(!1) and ~M0(!1) are examples

of low-pass lters (integration lters), whereasM1(!1) and ~M1(!1) are examples

of high-pass lters (dierentiation lters). So, we may interpret a two dimensional lter in our 2D lter bank, as a cascade of integration or dierentiation in thex1

-direction followed by integration or dierentiation in the x2-direction. However,

it seems to be unnatural to apply tensor product wavelets in problems where the Laplacian is involved. With respect to the Laplacian, scaling functions and wavelets which are radial symmetric should be preferred. But, it can be shown that in the lter bank setting these scaling functions and wavelets do not exist. A compromise could be the scaling functions and wavelets introduced by Cohen and Schenkler in (2]). Cohen and Schenkler consider scaling functions and wavelets on an hexagonal grid, T say (cf. Figure 3), which is generated by the vectors

e1 = (10),e2= ( ; 1 2 p 3 2 ). ande 3 = ;e 1 ;e 2 = ( ; 1 2 ; p 3 2 ) (cf Figure 4).

The spaceV0consists of all those continuous functions inL 2(IR

2), from which

the restrictions to the triangles of T coincide with linear functions. The

under-lying scaling function (x) is the function, which vanishes outside the hexagon

H (cf Figure 4). This function is also known as the Courant-Hilbert function.

Evidently, the translations are also adapted to the hexagonal grid, so we de ne

nj(x1x2) = 2 j(2jx 1 ;n 1e12 jx 2 ;n 2e2) (n= (n1n2) 2ZZ 2):

(15)

x1

x2

Figure 3: Hexagonaal grid T

e1

e2

e3

Figure 4: The hexagonH

An elegant representation of the scaling function  can be given in terms of convolution integrals. In order to obtain such a representation, let U be the

parallelogram U =fx=e 1+e2  0<  <10<  <1 g

and let U(x) be the characteristic function of U, i.e.

U(x) =

(

1 (x2U)

0 (x =2U):

Then, (x) may be represented by means of the following convolution type integral: (x) =Z 1 0 U(x;te 3)dt: (25)

(16)

^ (!) = 1;e ;i!e 1 i!e 1 1;e ;i!e 2 i!e 2 1;e ;i!e 3 i!e 3 (! 2IR 2):

>From this, it easily follows that ^ (!) = 1 8(1 +e ;i!e 1= 2)(1 +e;i!e 2= 2)(1 +e;i!e 3= 2) ^(!=2): Hence (cf (15)) M0(!) = 1 8(1 +e ;i!e 1)(1 +e ;i!e 2)(1 +e ;i!e 3): (26)

The scaling function is invariant with respect to rotation over an angle of 2=3, i.e.

(Rx) =(x) (x2IR 2)

where R is the rotation operator given by the matrix

R=  1 2 ; 1 2 p 3 1 2 p 3 1 2 ! :

Cohen and Schrenkler ( cf Cohen]) succeeded to compute wavelets l

cor-responding to , which are global rotational invariant with respect to R in the sense that

2(x) = 1(Rx) 3(x) = 2(Rx) (x 2IR

2): (27)

They presented dierent examples for such functions 1(x). In this chapter

we will use the function 1(x) which is associated with the following function

M1(!). M1(!) = ; 12 + 6e;i!e 1 ;24ei! e 1 + 4e 2i!e 1 + 2e 3i!e 1 + 5e ;i!e 2 + 5ei!e 2 ;2e ;2i!e 2 ;2e 2i!e 1 ;e ;3i!e 2 ;e 3i!e 2 ;24e i!(e 1 +e 2 )+ 4e2i!(e 1 +e 2 )+ 2e3i!(e1+e 2 )+ 6e;i!(e 1 +e2) ;e ;i!(e1+2e2) ; e;i!(e1+3e2) ;e ;i!(e1;e2) ;e ;i!(e1;2e2)+ei!(e1;e 2 ) ;ei! (e1;2e 2 )+ 2ei!(2e 1 ;e 2 )+ei!(e 1 +2e 2 ) ;ei! (e 1 +3e 2 )+ 2ei!(2e 1 +3e 2 )+ 2ei!(3e 1 +2e 2 )+ 2ei!(3e 1 +e 2 )+ 4ei!(2e 1 +e 2 )  =64: (28)

Because of (27), the functionM2(!) is obtained fromM1(!) by replacinge1

by e2,e2 by ;e

1 ;e

2. Subsequently, the functionM3(!) is obtained fromM1(!)

by replacing e1 by ;e

1 ;e

(17)

part of the lterbank are known for this moment. For the decomposition part, we need the duals ~Ml(!)(l = 0123). Again, Cohen and Schenkler presented

in Cohen] the following formula: 32 ~M0(!) = 14 + 6e ;i!e 1 + 6ei! e 1 ;e ;2i!e 1 ;e 2i!e 1 + 6e ;i!e 2 + 6ei!e 2 ;e ;2i!e 2 ;e 2i!e 2 + 6e ;i!(e 1 +e 2 )+ 6ei!(e 1 +e 2 ) ; e;2i!(e 1 +e 2 ) ;e 2i!(e 1 +e 2 ) ;2e ;i!(e 1 ;e 2 ) ;2ei! (e 1 ;e 2 ) ;2e ;i!(e 1 +2e 2 ) ;2e i!(e 1 +2e 2 ) ;2e ;i!(2e 1 +e 2 ) ;2e i!(2e 1 +e 2 ) (29) Finally, the decomposition lter ~M1(!) is given by

4 ~M1(!) = 1 ;ei! e 1 ;ei! (e 1 +e 2 )+ei!(2e 1 +e2) (30)

Evidently, the functions ~Ml(!)(l = 23) can be derived from ~M1(!) by

re-placing the vectors ej in an appropriate way.

Now, our lter bank is completed, since the lter coecients are known.

Example

As an application we consider the smooth Gaussian functionf(x1x2) =e ;(x 2 1 +x 2 2 ),

which is sampled at the points

h(n1e1+n2e2), (n1n2 =

;N;N+ 1::: N),

withN = 20 andh= 0:05. In fact the spaceV0now consists of the continuous

piecewise linear functions on the hexagonal grid hT. However, sampling the

function f(x) at the grid hT coincides with sampling the functionf(hx) at the

gridT. So, we may apply our lter bank to the matrixan

0 =f(h(n1e1+n2e2)).

At the grid points of hT, where a sampled value is not available, it is assumed

that an0 = 0. Therefore, we have extended an0 to the whole grid by zero

padding. This zero padding will introduce boundary errors in the reconstructed values as will be clear by the next tables. The input of the decompostion part of the lter bank is a 4141 matrix of an

0 values. The output consists of four

2020 matrices representing the values of an ;1,d 1 n;1,d 2 n;1,d 3 n;1respectively.

In the four output matrices, we replace all the elements for which the absolute value is less then a certain treshold value by zeros. Then the reconstruction part of the lter bank is applied. For an exact reconstruction, the output of the reconstruction part (the reconstruction matrix) should be equal to the input of the decomposition part. Because of tresholding and zeropadding this will not be the case. In the third column of the tables below we have plot the error"1, which

is the maximum absolute error in the reconstruction matrix, compared to the original matrix. As will be clear from the tables, this error is mainly caused by

(18)

zeropadding. The error due to zeropadding will occur at the "boundary" of the matrix. Therefore, the fourth column contains the absolute error ("2), without

taking into account the rst three colums and rows and the last three columns and rows of the reconstruction matrix. In "2, the error due to zeropadding is

eliminated, so this error is mainly caused by tresholding. It can be shown that it is comparable to the treshold value. In the rst column the dierent applied treshold values are listed. The second column gives the ratio of the total number of nonzero elements of the four output matrices (after tresholding) and the size of the input matrix (an0. This ratio is the compression factor.

treshold compr. factor "1 "2

0.01 0.26 3.8 10;1 4.5 10;3 0.001 0.42 3.8 10;1 3.2 10;3 0.0001 0.90 3.8 10;1 1.8 10;4 0 0.95 3.8 10;1 0 f(x1x2) = exp( ;x 2 1+x 2 2), h= 0:05,N = 20

treshold compr. factor "1 "2

0.01 0.32 8.7 10;1 6.9 10;3 0.001 0.77 8.7 10;1 2.2 10;3 0.0001 0.91 8.7 10;1 0 0 0.95 8.7 10;1 0 f(x1x2) = sin( jx 1 ;0:1j+jx 2+ 0:1 j), h= 0:05,N = 20

It is clear that small treshold values have negative in#uence on the compres-sion factors. Since the approximation error by linear interpolation has order h2

and the error due to tresholding is comparable to the treshold value, to use a treshold value of order h2 seems to be a good choice. In the next table, we used

the function f(x1x2) = e ;x 2 1 ;x 2

2 again, but h= 0:05 is replaced byh = 0:1 and

N = 20 is replaced byN = 10.

treshold compr. factor "1 "2

0.01 0.27 3.8 10;1 1.9 10;2 0.001 0.80 3.8 10;1 2.3 10;3 0.0001 0.91 3.8 10;1 4.7 10;5 0 091 3.8 10;1 0 f(x1x2) = exp( ;x 2 1+x 2 2), h= 0:1, N = 10

(19)

4 The computation of a particular solution

In this section, we will focus on the problem, how to compute a particular solution of the dierential equation u(x) =f(x) (x2D), for a given functionf(x) and

a bounded domain D. First, the function f will be approximated by piecewise

linear interpolation on a part of the hexagonal gridhT, which contains the given

domain D. Let f

0(x) be the approximating function, then

f(x) f 0(x) = X n1n2 f(hn1e1+hn2e2)(x=h ;n 1e1 ;n 2e2)

where(x) is our scaling function for the hexagonal wavelets. Due to linear interpolation, the approximation error onDis of orderh

2, under

the assumption that f is a continuous function. In order to apply compression, as demonstrated in the previous section, the decomposition part of the lter bank (cf Figure 1) associated to the hexagonal wavelets must be used. The four output matrices of the decomposition part contain the coecients for the decompostion of f0 as given by (22). In this decomposition the function f0 is

written as a linear combination of dilated and translated versions of the function

(x), and the three mother wavelets 1(x), 2(x) and 3(x). Tresholding the

coecients in this combination which are small compared to a given treshold value will introduce a second error, which can be controlled by computing the reconstruction error. This error is comparable to the treshold value. Because of (12), to nd a particular solution corresponding to the wavelets, it is sucient to nd a particular solution correponding to functions of the type(2

;l

h x;n 1e1

;

n2e2). Ifp(x) is a particular solution corresponding to the scaling function(x),

then 4lh2p(2;lx=h ;n

1e1 ;n

2e2) is a particular solution for(2 ;lx=h

;n 1e1

;

n2e2). Therefore we only have to computep(x).

Evidently, a particular solution p(x) can be expressed by the convolution integral: p(x) =ZZ IR 2 ()G(x;)d: (31) where G(x) = 1 2 log

jxj is the fundamental solution correponding to the

Laplacian.

Our main task is to evaluate the integral (31). By substituting expression (25) into (31), we get

(20)

p(x) =Z 1 0 ZZ IR 2 U(;te 3)G(x ;)d dt= Z 1 0 ZZ IR 2 U()G(x;;te 3)d dt =Z 1 0 ZZ U G(x;;te 3)d dt: Now, we replace  = (12) by  0 = (0 1 0 2) such that  =  0 1e 1 + 0 2e 2. Hence, 1= 0 1 ; 1 2 0 2 and  2 = 1 2 p 30 2: By setting 0

3=t, we arrive at the representation

p(x) = 1 2 p 3R 1 0 R 1 0 R 1 0 G(x ; 0 1e 1 ; 0 2e 2 ; 0 3e 3)d 0 1d 0 2d 0 3: (32)

The next natural step is to represent x 2 IR

2 in terms of e

1 e2 and e3 as

follows:

x=1e1+2e2+3e3

1+2+3= 1:

The numbers 1, 2 and 3 are called barycentric coordinates of x with

respect to e1,e2 and e3. By substituting 0 1 =  1 ;u 0 2 =  2 ;v and  0 3 =  3 ;w, the following relation is obtained. p(x) = p 3 8 Z  3 3 ;1 Z  2 2 ;1 Z  1 1 ;1 log(jue 1+v e2+we3 j 2)dudv dw: (33) Since jue 1 +v e2+we3 j 2 = 1 2((u ;v) 2+ (u ;w) 2 + (v ;w) 2), a function

F(uvw) for which

@ @w @v@ @u F@ (uvw) = log((u;v) 2+ (u ;w) 2+ (v ;w) 2)

will deliver our p(x) =p(1e1+2e2+3e3) as follows

p(x) = p 3 8 (;log 2 + (I;E 3)(I ;E 2)(I ;E 1)F( 123)): (34) Here, I is the identity operator andEa denotes the backward shift operator

with respect to the variabele a, for instanceE1F(

123) =F(1 ;1

(21)

In this way, our function p(x) is expressed as a third order nite dierence of a function F(uvw) on a cube in IR3. A function F(uvw), satisfying the

previous partial dierential equation can be found by repeated integration. An explicit formula is given by the expression

F(uvw) =h1(uvw) + h2(uvw)(log((u ;v) 2+ (u ;w) 2+ (v ;w) 2) ; 11 3) (35) where, h1(uvw) = p 3 6 (u;v) 3 arctan(u+v;2w p 3 (u;v) ) + p 3 6 (u;w) 3 arctan(u+w;2v p 3 (u;w)) + p 3 6 (v;w) 3 arctan(v+w;2u p 3 (v;w))  (36) and h2(uvw) = 1 12(2v ;u;w)(2u;v;w)(2w;u;v): (37)

Our conclusion is that p(x) can be found by doing elementary computations. In BEM also normal derivatives at some boundary points of ; of are needed. Formulas for px1(x) =

@p

@x1(x) and px 2(x) =

@p

@x2(x) are given below. In these

formulas, we use again the barycentric coordinates1,2and3ofxwith respect

to e1,e2 and e3 and functions Fk(123), (k = 12), which play a similar role

as the functionF(123) forp(x). These functions are given by the formulas:

F1(uvw) = ; 2p 3(u;v) 2 arctan(u+v ;2w p 3(u;v) ) + 2p 3(u;w) 2 arctan(u+w ;2v p 3(u;w) ) + ((u;v) 2+ (u ;w) 2 ;2(v;w) 2) log(((u ;v) 2+ (u ;w) 2+ (v ;w) 2))  =4: F2(uvw) = ; 1 2 (u;v) 2 arctan(u+v ;2w p 3(u;v) ) + 1 2 (u;w) 2 arctan(u+w ;2u p 3(u;w) ) + ( v;w) 2 arctan(v+w ;2u p 3(v;w) ) + p 3 4 (v;w)(2u;v;w)(3;log((u;v) 2+ (u ;w) 2+ (v ;w) 2)):

(22)

Finally, the partial derivatives ofp(x) are given by px1(x) = p 3 8 (I;E 1)(I ;E 2)(I ;E 3)F 1(123) px2(x) = p 3 8 (I;E 1)(I ;E 2)(I ;E 3)F 2(123):

Here 1, 2 and 3 are the barycentric coordinates of x with respect to e1,

e2 ande3.

Example

As a numerical test we solve the Dirichlet problem for

f(x1x2) = 14 (x 2 1+x 2 2 ;1)e ;x 2 1 ;x 2 2 g(x1x2) = e ;1

on the unit circle x2 1 +x 2 2  1 with boundary ; : x 2 1+x 2 2 = 1. The exact solution equals u(x1x2) = e ;x 2 1 ;x 2

2. We solve this problem by using a BEM

method having 10 quadratic boundary elements and used the hexagonal wavelets for deriving a particular solution by the method described in this section with meshsize h = 0:05, n = 50 and treshold value equals 0.001. The compression factor is equal to 0.25 .The dierence between the exact solution and the BEM solution at certain points in the circle are tabulated below

x1 x2 wavelet exact

0.0 0.0 0.165947 0.165979 0.25 0.0 0.198489 0.198554 0.75 0.0 0.22557 0.225667 0.90 0.0 0.00000 0.000000

We conclude that the approximation is of order 10;4, which is smaller then

h2. However, in general we may not expect to have an approximation orderhp,

with pless then 2.

References

1] C.K. Chui, An introduction to wavelets, Academic Press, Inc (1992)

2] A. Cohen, J.M. Schenkler (1993): Compactly supported bidimensional wavelet bases with hexagonal symmetry. Constr. Approx. (1993) 9: 209-236 3] I. Daubechies, Ten lectures on wavelets, SIAM (1991)

(23)

4] M.A. Golberg, The evaluation of particular solution in the bem, Boundary elements communications 6 (1995), 99-106

5] P.W. Partridge, C.A. Brebbia, L.C. Wrobel, The dual boundary element method, Computational Mechanics Publications, Elsevier (1992)

6] G. Schmidt, H. Strese, BEM for Poisson equation, Engineering analysis with boundary elements 10 (1992), 119-123

7] C. Sidney Burrus, Ramesh A. Gopinath, Haito Guo, Wavelets and Wavelets Transforms, Prentice Hall (1998)

8] G. Strang, T.Nguyen, Wavelets and Filter Banks,Wellesly-Cambridge Press (1996)

Referenties

GERELATEERDE DOCUMENTEN

If the research can describe the influence of HIV/AIDS, the role that the global and the local context plays in it and if the research can help the

The Constitution grants the right for the various political parties to appoint and dismiss members of the parliament and therefore, if the President acts as the leader of the

4 Godsbegrip  4.1 Doel van die hoofstuk 

This research seeks to establish the political role that the City Press defined for its black journalists in post-apartheid South Africa, and the role played by

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers) Please check the document version of this publication:.. • A submitted manuscript is

• The final published version features the final layout of the paper including the volume, issue and page numbers.. Link

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers) Please check the document version of this publication:.. • A submitted manuscript is

The condition number of the matrices A (circles) and G (squares), corre- sponding to the Laplace equation with mixed boundary conditions and Dirichlet boundary conditions