A meshless polyharmonic-type boundary interpolation method for solving boundary integral equations

Engineering Analysis with Boundary Elements 28 (2004) 1207–1216 www.elsevier.com/locate/enganabound

A meshless polyharmonic-type boundary interpolation method for solving boundary integral equations Csaba Ga´spa´r Department of Mathematics, Sze´chenyi Istva´n University, P.O. Box 701, H-9007 Gyo¨r, Hungary Received 15 July 2002; revised 5 March 2003; accepted 20 April 2003 Available online 10 March 2004

Abstract A boundary interpolation technique is introduced based on multi-elliptic partial differential equations. The interpolation problem is converted to a special higher order partial differential equation which is completely independent of the geometry of the original problem. Based on this interpolation method, meshless methods are constructed for the 2D Laplace – Poisson equation. The presented approach makes it possible to avoid solving large and dense interpolation equations. The auxiliary higher order partial differential equation is solved by robust, quadtree-based multi-level methods. The results can be easily generalized to 3D problems as well. q 2004 Elsevier Ltd. All rights reserved. Keywords: Meshless methods; Scattered data interpolation; Biharmonic interpolation; Polyharmonic interpolation; Quadtree

1. Introduction The success of the boundary element method is mainly based on two important features. First, the boundary element method reduces the dimension of the original problem by one, thus, the number of the introduced unknowns is much less than that of the traditional domain type methods. Second, the domain mesh generation is often avoided, which is generally the most difficult task when solving a partial differential equation. However, to solve the associated boundary integral equations, a boundary mesh is still needed in general. On the other hand, one must pay a certain price for these advantageous features. The standard discretization techniques of the boundary element method (point collocation, least squares, Galerkin method, etc.) lead to algebraic equations with fully populated and non-selfadjoint matrices, while the classical domain type methods produce sparse matrices. This may result in numerical difficulties: a traditional Gaussian elimination requires OðN 3 Þ algebraic operations (where N denotes the number of boundary nodes). If N is large, this may be inadmissibly high and the solution procedure may become unstable. Moreover, if the original partial differential equation is inhomogeneous and/or nonlinear, domain integrals also appear. E-mail address: [email protected] (C. Ga´spa´r). 0955-7997/$ - see front matter q 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2003.04.001

The evaluation of the domain integrals should be performed without using any domain mesh: otherwise, one of the main advantages of the boundary element method would be lost. These disadvantages are well known and a number of sophisticated techniques have been developed to overcome the difficulties indicated above. The computational cost can be significantly reduced by applying special iterative and multi-grid techniques [10]. Another way to reduce the computational cost is to apply the fast multipole method of Rokhlin [25] in the evaluation of the appearing boundary integrals. This method can be combined with multi-level techniques [12]. It is also possible to apply a domain decomposition technique, especially when it is used with boundary-type preconditioners [20,21]. In addition to this, domain decomposition techniques make the solution procedure easily parallelizable. As to the treatment of the domain integrals, the most popular technique is the dual reciprocity method of Partridge and Brebbia [24], which converts the appearing domain integrals to boundary integrals. This approach requires a set of points inside the domain and is based on a scattered data interpolation. The interpolation points need not have any structure, so that domain mesh generation is still avoided. (This is not the case on the boundary.) Depending on the choice of the interpolation technique, a series of such methods can be derived. It should be emphasized, however, that in a significant part of these

1208

C. Ga´spa´r / Engineering Analysis with Boundary Elements 28 (2004) 1207–1216

methods (when globally supported basis functions are used), the interpolation subproblem suffers from the same computational disadvantages than the Boundary Element Method itself as pointed out earlier. It leads to a large and often non-selfadjoint algebraic system with a fully populated matrix which is inherently ill-conditioned when the interpolation point set contains points which are located too near to each others. These difficulties can be significantly reduced if the biharmonic interpolation is applied [9,13]. In this approach, the interpolation problem is converted to the biharmonic equation which is supplied with the interpolation equations as special boundary conditions. Pointwise boundary conditions (at discrete points) make the usual second-order partial differential equations ill-posed, but this is not the case at higher order partial differential equations as the biharmonic equation. The method is closely related with the augmented thin plate splines, but no large and illconditioned system has to be solved. Instead, it is the biharmonic interpolation problem that should be solved, which is, however, completely independent of the original domain V and also its boundary. Using robust, multi-level quadtree (QT)-based solvers, the computational cost of the interpolations subproblem can be reduced to OðM log MÞ, where M denotes the number of the interpolation points scattered in V: For details, see Refs. [13,14]. Even if some of the above outlined methods is used, a boundary mesh is still necessary in order to evaluate the appearing boundary integrals. This task is easy in 2D problems but not in 3D. To completely get rid of the mesh generation, several techniques have been developed. The main idea is again the application of the scattered data interpolation. As a special case, the idea of the indirect BEM has also returned into practice with the important modification that the point sources are located outside the domain V and not along the boundary (see Refs. [1,3]). Unfortunately, only very smooth solutions can be obtained in this way, and the source points should be located not too far from the boundary of V : otherwise, the discrete problem becomes ill-conditioned quite soon. The fast multipole summation idea can be applied again in order to reduce the computational complexity [2]. In this paper we apply an interpolation approach which is similar to that of [9,13]. The appearing domain integrals are evaluated using a biharmonic-type scattered data interpolation based on discrete points inside the domain. The homogeneous problem is solved by a boundary interpolation using another fourth-order partial differential operator. Thus, neither domain nor boundary mesh generation is necessary: in addition to this, the problems of large, full algebraic systems are also avoided. In Section 2, we summarize the main tools and results concerning the scattared data interpolation based on polyharmonic and other poly-elliptic partial differential operators including the questions of numerical realization in QT context. In Section 3, we detail the above indicated meshless method and present some simple numerical

examples. The construction will be presented through the 2D model problem: Du ¼ f

ð1Þ

in V

supplied with Dirichlet boundary condition: ul›V ¼ u0 : Here V , R2 denotes a bounded, sufficiently smooth domain. To discretize the problem (1), we need a set of distinct points scattered inside the domain: x1 ; …; xM [ V and another set of distinct points scattered along the boundary: xMþ1 ; …; xMþN [ ›V: The function f is given by the values f ðx1 Þ; …; f ðxM Þ: Note that, though the model problem is a 2D one, the results and the algorithms can easily be generalized to 3D problems as well, with natural modifications (replacing the appearing 2D fundamental solutions with 3D fundamental solutions; QTs with octtrees and so forth.)

2. Scattered data interpolation based on iterated elliptic equations 2.1. Scattered data interpolation Let V0 , R2 be a bounded, smooth domain (not necessarily identical to the domain V in Eq. (1)). Let x1 ; x2 ; …; xN [ V0 be distinct points and let f1 ; f2 ; …; fN [ R be given values. The problem is to find a (sufficiently regular) function f : V0 ! R such that it satisfies the interpolation conditions: f ðxk Þ ¼ fk

ðk ¼ 1; 2; …; NÞ

ð2Þ

Of course, one should specify a proper function space in which the solutions of the problem are sought. To investigate the properties of different interpolations, the above problem is reformulated as follows: Problem 1. Let X; Y be Banach spaces, X , CðV 0 Þ (the space of the functions continuous on the compact set V 0 supplied with the maximum norm). Suppose that X , Y with continuous injection. Find an interpolation operator A : X ! Y with respect to the set {x1 ; x2 ; …; xN }; i.e. a bounded linear operator such that ðAf0 Þðxk Þ ¼ f0 ðxk Þ

ðk ¼ 1; 2; …; NÞ

holds for every function f0 [ X: In constructing interpolation operators, the goal is that A should approximate the identity operator I as accurately as possible. The interpolation is said to be of order s $ 0 (with respect to the norms of X and Y), if kAf0 2 f0 kY # C0 hs kf0 kX holds for every f0 [ X and every finite point set {x1 ; x2 ; …; xN } , V0 ; where the constant C0 is independent of both f0 and the interpolation point set. Here h denotes

C. Ga´spa´r / Engineering Analysis with Boundary Elements 28 (2004) 1207–1216

the separation distance: h U sup min kx 2 xk k

ð3Þ

x[V0 1#k#N

2.2. Interpolation by radial basis functions Perhaps the most popular method for the scattered data interpolation problem is the method of radial basis functions (RBFs, see, e.g. Franke [6], Franke and Schaback [7], Golberg and Chen [16,17]). in which the interpolation function f U Af0 is defined in the form: N X f ðxÞ ¼ aj Fj ðx 2 xj Þ ð4Þ j¼1

where F1 ; F2 ; …; FN are predefined radial basis functions (circularly symmetric functions), i.e. the values of Fj ðxÞ depend on the norm r ¼ kxk only. The a priori unknown coefficients aj are determined by solving the interpolation equations N X aj Fj ðxk 2 xj Þ ¼ f0 ðxk Þ ðk ¼ 1; 2; …; NÞ ð5Þ j¼1

Depending on the choice of the applied RBFs, various methods can be constructed, e.g. qffiffiffiffiffiffiffiffiffi † Fj ðrÞ U c2j þ r 2 (Hardy’s method of multi-quadrics, with scaling parameters cj ; † Fj ðrÞ U r 2 log r (thin plate splines), or, more generally: † Fj ðrÞ U r 2m log r (polyharmonic splines) with a positive integer parameter m; † Fj ðrÞ U 1 þ r (applied often in dual reciprocity method, see Partridge and Brebbia [24]); † Fj ðrÞ U e2cj r (Gaussian functions) and so forth. The above RBFs are globally supported. Consequently, all of these methods lead to a large and full system of equations of type Eq. (5). These systems are inherently illconditioned if the distance of some interpolation points is too small (this results in nearly identical rows in the matrix). This is a serious drawback of these methods, though they may exhibit excellent approximation properties. It is known for example, that the method of multiquadrics is exponentially convergent with respect to the maximum norm provided that f0 is smooth enough [18,22]. A possible way to avoid full algebraic systems is the use of compactly supported RBFs as proposed by Wendland [26] and/or to apply multi-level solvers [4,5]. We will follow, however, another way, which is similar to the thin plate splines but without explicitly using these RBFs. 2.3. Interpolation methods based on higher order pdes The thin plate spline interpolation functions (as well as the augmented thin plate spline interpolation completed by linear polynomials) solve the biharmonic equation

1209

everywhere in V0 ; except the interpolation points. Similar observation can be made for the piecewise linear interpolation and also for the cubic splines in 1D: they satisfy the ordinary differential equation d2 f =dx2 ¼ 0 and d4 f =dx4 ¼ 0, respectively (except at the interpolation points). To generalize the idea, let L be a (carefully chosen) partial differential operator, and define the interpolation function f U Af0 as a solution of the following pde: Lf ¼ 0

in V0 \{x1 ; x2 ; …; xN }

ð6Þ

supplied with the interpolation conditions as special boundary conditions: f ðxk Þ ¼ f0 ðxk Þ

ðk ¼ 1; 2; …; NÞ

ð7Þ

Along the boundary ›V0 ; a usual regular boundary condition, e.g. Dirichlet boundary condition can be imposed. It should be pointed out that the boundary condition taken along ›V0 plays a minor role in the interpolation: the crucial point is the proper choice of the differential operator L: Since pointwise boundary conditions at discrete points are not quite usual in practice, a careful analysis must be carried out. The fundamental solution of the operator L will play an important role in the following. A few applicable operators together with the corresponding fundamental solution are listed in Table 1. The parameter c is a positive scaling constant, K0 , K1 denote the usual Bessel functions of the third kind. Note that the fundamental solutions of D and D 2 c2 I have a logarithmic singularity at the origin. The other fundamental solutions are continuous everywhere. Note also that, in the last two items, the singularities cancel out. In the next subsections we briefly summarize the main results of this kind of interpolation. Details as well as proofs can be found in Refs. [13 – 15]. 2.3.1. Harmonic interpolation The interpolation of the form (6) – (7) fails if L is a familiar second-order operator, e.g. the Laplacian. To illustrate this phenomenon, consider the following example. Let f0 ðx; yÞ U sin2 px sin2 py

ð8Þ

and seek the function f as a solution of the Laplace equation: Df ¼ 0 in V0 \{x1 ; x2 ; …; xN } supplied with the boundary Table 1 2D partial differential operators and their fundamental solutions Operator

Fundamental solution (FðrÞÞ

D DD D3 ðD 2 c2 IÞ ðD 2 c2 IÞ2 DðD 2 c2 IÞ DðD 2 c2 IÞ2

1 2p log r 1 2 8p r log r 4 1 128p r log r 1 2 2p K0 ðcrÞ r 4pc K1 ðcrÞ 1 2 2pc 2 ðK0 ðcrÞ þ log crÞ 1 r ðK0 ðcrÞ þ log crÞ þ 4pc 3 2pc4

K1 ðcrÞ

C. Ga´spa´r / Engineering Analysis with Boundary Elements 28 (2004) 1207–1216

1210

Direct interpolation problem: Find a function u [ W such that ð11Þ DDðu þ f0 Þ ¼ 0 in V0 \{x1 ; x2 ; …; xN } is valid (in the sense of distributions). Variational problem. Find a function u [ W such that for every w [ W : ð Dðu þ f0 ÞDw dV ¼ 0 ð12Þ V0

Fig. 1. Harmonic interpolation with logarithmic singularities.

condition f l›V0 ¼ 0 and the interpolation conditions f ðxk Þ ¼ f0 ðxk Þ ðk ¼ 1; 2; …; NÞ: In this example, V0 is the unit square and the interpolation points are quasi-randomly scattered in V0 with N ¼ 100: The Laplace equation is discretized by an equidistant 64 £ 64 grid. The discrete solution is shown in Fig. 1. One can easily observe the (logarithmic-type) singularities at the interpolation points, which destroys the approximation of the interpolation (the relative error was 13.49% with respect to the discrete L2 -norm using the above discretization). If the grid size tends to zero, the singularities become larger and larger and the corresponding discrete interpolation functions do not converge. The reason of the appearance of the singularities is that the fundamental solution of the Laplacian is not continuous at the origin, therefore the above harmonic interpolation problem cannot be a well-posed problem in any function space contained in CðV 0 Þ provided that f0 is not identically zero at the interpolation points. Note that, in 1D, the corresponding differential operator d2 =dx2 does have a continuous fundamental solution ð1=2lxlÞ; and the 1D harmonic interpolation results in the familiar piecewise linear interpolation. In 2D and 3D, however, higher order partial differential operators are needed. 2.3.2. Biharmonic interpolation Defining the operator L by L U DD; we obtain the biharmonic interpolation problem: DDf ¼ 0 in V0 \{x1 ; x2 ; …; xN } f l›V0 ¼ 0;

ð9Þ

›f l ¼ 0 ðboundary conditionsÞ ›n ›V0

f ðxk Þ ¼ f0 ðxk Þ ðk ¼ 1; 2; …; NÞ ðinterpolation conditionsÞ To analyze the solvability and the approximation properties, we give a variational form of this problem. Let us introduce the following subspace of the Sobolev space H02 ðV0 Þ : W U {u [ H02 ðV0 Þ : uðx1 Þ ¼ … ¼ uðxN Þ ¼ 0}

ð10Þ

Due to the well-known imbedding theorems, this subspace is closed in H02 ðV0 Þ: Now we reformulate the biharmonic interpolation problem (9) as follows:

is satisfied. It can be shown (see Refs. [14,15]) that the direct and the variational problems Eqs. (11) and (12) are equivalent and have a unique solution for arbitrary function f0 [ H02 ðV0 Þ; i.e. the interpolation conditions (which can be considered as unusual boundary conditions taken at discrete points) do not destroy the well-posedness of the boundary value problem of the biharmonic equation in Sobolev space H02 ðV0 Þ: Having solved either Eqs. (11) or (12), the biharmonic interpolation function is defined by f U u þ f0 : The biharmonic interpolation is closely related with the augmented thin plate splines. More precisely, the biharmonic interpolation function f defined by the solutions of Eqs. (11) and (12) is uniquely represented in the following form: N X f ðxÞ ¼ wðxÞ þ bj Fðx 2 xj Þ ð13Þ j¼1

where w is a function which is biharmonic everywhere in V0 and F is the biharmonic fundamental solution: FðrÞ ¼ 1 2 8p r log r: This shows that the interpolation function can be expressed as a sum of a regular function w and a linear combination of weakly singular radial basis functions. In the case of the augmented thin plate splines, the function w is linear. It should be pointed out that once the interpolation function f has been determined, the coefficients b1 ; b2 ; …; bN and the function w can also be determined without solving any additional system of equations. Indeed, the representation (13) immediately implies that (in the sense of distributions): N X DDf ¼ b j d xj ð14Þ j¼1

where dxj denotes the Dirac distribution concentrated to the point xj : From a proper discretization of Eq. (14), the coefficients bj can be directly computed. Roughly speaking, bj is the ‘intensity’ of the weak singularity of f at the point xj : As to the approximation of the biharmonic interpolation, it has been proved that the operator A defined by Af0 U f is second-order with respect to the L2 ðV0 Þ-norm and firstorder with respect to the H 1 ðV0 Þ-norm, i.e. the estimations kAf0 2 f0 kL2 ðV0 Þ # C1 h2 kf0 kH02 ðV0 Þ ; kAf0 2 f0 kH 1 ðV0 Þ # C2 hkf0 kH02 ðV0 Þ hold for every function f0 [ H02 ðV0 Þ; where h denotes the separation distance defined by Eq. (3), and the constants C1 ;

C. Ga´spa´r / Engineering Analysis with Boundary Elements 28 (2004) 1207–1216

Fig. 2. Biharmonic interpolation: no singularities occur.

C2 are independent of both f0 and the interpolation point set {x1 ; x2 ; …; xN }: As an illustrative example, consider again the function (8) of the previous example. The set of interpolation points remains the same but a biharmonic interpolation is applied instead of the harmonic interpolation. The result is plotted in Fig. 2. The interpolation became much better, the relative error (measured in the discrete L2 -norm) was reduced to 1.43%. 2.3.3. Polyharmonic and other poly-elliptic interpolations A natural generalization of the biharmonic interpolation is the use of higher powers of the Laplacian, i.e. polyharmonic operators. Let m . 2 be an integer and L U Dm ; then we arrive at the m-harmonic interpolation problem:

1211

which is closed in H0m ðV0 Þ; if m $ 4; and so forth. It is also possible to use the subspace W0;m > W1;m ; which corresponds with the interpolation condition when both f and Df are prescribed at the interpolation points. A further possibility is to prescribe the values of f at certain points and the values of Df at the remaining points. This strategy fits well a natural meshless technique for the Poisson equation as will be shown in Section 3. In all these cases, analogous solvability and representation theorems can be derived with the important modification that the representation theorems may contain also the derivatives of the fundamental solution of the applied differential operator. The order of approximation grows with the value of m : in practice, however, we have found that m should not exceed the value 4. Otherwise, the numerical implementation becomes too expensive. Another way to generalize the biharmonic interpolation idea is to replace the Laplace operator with the negative definite elliptic operator D 2 c2 I; where c . 0 is a scaling constant. The simplest case of this approach is the use of the bi-Helmholtz operator L U ðD 2 c2 IÞ2 : The corresponding theorems are quite similar to those of the biharmonic interpolation. The representation formula is based on the fundamental solution of the operator ðD 2 c2 IÞ2 ; i.e.

FðrÞ ¼

r K ðcrÞ 4pc 1

ð16Þ

If the values of Df are to be prescribed at the interpolation points, the corresponding subspace is:

The approximation properties are also similar to those of those biharmonic operator. The essential difference between the biharmonic and bi-Helmholtz interpolation is that while 1 2 the biharmonic fundamental solution FðrÞ ¼ 8p r log r is a globally supported function and lFðrÞl ! þ1 when r ! þ1; the latter fundamental solution (16) is a rapidly decreasing function: lFðrÞl ! 0 quickly when r ! þ1: Though theoretically it is also globally supported, from computational point of view it can be considered as if it were compactly supported. The size of the ‘computational support’ can be controlled by the parameter c: As an important consequence of this feature is that the interpolation works in the case when V0 ¼ R2 and PN still p f ðxÞ ¼ j¼1 bj Fðx 2 xj Þ (the representation formula does not contain the regular function w). This expression is exactly a classical RBF-type interpolation formula based on the fundamental solution of the bi-Helmholtz operator. Consequently, the above theory immediately results in existence, uniqueness and approximation theorems to that type of interpolation as well. Setting the constant c large enough makes the matrix of the interpolation equations sparse. However, if c is chosen too large, the computational support of F becomes smaller than the separation distance of the interpolation points, which generates singularities at the interpolation points, and the interpolation breaks down. Finally, consider a mixed Laplace – Helmholtz-type interpolation. Here the operator L is defined by

W1;m U {u [ H0m ðV0 Þ : Duðx1 Þ ¼ Duðx2 Þ… ¼ DuðxN Þ ¼ 0}

L U DðD 2 c2 IÞ

Dm f ¼ 0 in V0 \{x1 ; x2 ; …; xN } ðdifferential equationÞ

›l f l›V ¼ 0 ðl ¼ 0; 1; …; m 2 1Þ ðboundary conditionsÞ ›nl 0 f ðxk Þ ¼ f0 ðxk Þ ðk ¼ 1; 2; …; NÞ ðinterpolation conditionsÞ ð15Þ The boundary conditions can be replaced again with other regular boundary conditions. As to the interpolation conditions, now we have much greater degree of freedom. Since we will seek the solution of Eq. (15) in a subspace of H0m ðV0 Þ; due to the imbedding theorem, the solution of the problem (15) (if exists), not only continuous in V0 but all derivatives up to the ðm 2 2Þth order are also continuous. Therefore it is anticipated that not only the values but also the derivatives of f (up to the ðm 2 2Þth order) can also be prescribed at the interpolation points (Hermite-type interpolation). According to the choice of the interpolation conditions and the associated closed subspace of H0m ðV0 Þ; various direct and variational problems can be defined. For instance, the interpolation conditions in Eq. (15) correspond with the subspace W0;m U {u [ H0m ðV0 Þ : uðx1 Þ ¼ uðx2 Þ… ¼ uðxN Þ ¼ 0}

ð17Þ

1212

C. Ga´spa´r / Engineering Analysis with Boundary Elements 28 (2004) 1207–1216

where c . 0 is again a scaling parameter. Using the same closed subspace W , H02 ðV0 Þ as in the biharmonic case, similar theorems can be deduced. In the corresponding representation, however, the fundamental solution of the operator DðD 2 c2 IÞ appears:

FðrÞ ¼ 2

1 ðK0 ðcrÞ þ logðcrÞÞ 2pc2

ð18Þ

Note that F is continuous also at the origin, since the singularities of K0 and the logarithm function cancel out. Since the function K0 rapidly decreases, the function F becomes (approximately) equal to the harmonic fundamental solution far from the origin. Therefore this type of interpolation can be considered an ‘improved harmonic’ or ‘quasi-harmonic’ interpolation which does not generate singularities. This property makes it possible to construct a boundary meshless method for the Laplace equation as will be shown in Section 3. 2.4. Numerical implementation using quadtrees In the previous subsections, the scattered data interpolation problem was converted to an auxiliary fourth (or even higher) order partial differential equation. Since the scattered data interpolation is used to construct (meshless) solution methods for other partial differential equations, we have arrived at an apparent contradiction. In order to solve, e.g. the simple model problem (1), another and even more complicated partial differential equation, e.g. the biharmonic equation has to be solved. It seems that the computational complexity of the solution of the more complicated biharmonic equation will necessarily exceed that of the simpler model Poisson equation. However, there is an essential difference between the model problem and the associated interpolation problem which is converted to a higher order partial differential equation. Namely, the original domain V as well as its boundary ›V is completely missing from the latter problem. Instead, a larger domain V0 is needed, but it can be defined in a practically arbitrary way, so that it can be as simple as possible, e.g. a square or a circle. Note that, if the domain is simple, a number of well-known and highly efficient solvers are available (e.g. Fast Fourier Transform method or multigrid methods). What makes the solution computationally difficult is the possibly complicated geometry, and it is exactly this difficulty that can be circumvented by the presented approach. To solve the appearing biharmonic and similar interpolation problems, we used the multi-level QT-based method, which has proved a very efficient tool in quite different areas, e.g. numerical grid generation, databases, fast multi-pole method, computational flow modelling, etc. See, e.g. Greaves and Borthwick [19], Yiu et al. [27]. The QT algorithm is a recursively defined subdivision procedure of an initial square which results in a nonequidistant, non-uniform nested cell system as a computational grid. The structure of the cell system can finely be

controlled by a finite set of controlling points, so that local refinements can easily be generated. It should be pointed out that the algorithm is quite economic in the sense that it requires typically OðN log NÞ arithmetic operations only, where N denotes the number of the controlling points. In the interpolation methods described above, we used the simplest cell-centered box schemes to discretize the Laplace operator. (More sophisticated (and more exact) schemes can be defined by the method of moving least squares instead of box schemes (see, e.g. Morinishi [23])). Note that, it is sufficient to discretize the operator of the type ðD 2 c2 IÞ only, since the appearing equations can be split into such equations. For instance, the biharmonic equations DDf ¼ 0 is equivalent to the pair of Poisson equations Df ¼ g; Dg ¼ 0: Special care should be taken of the proper handling of the interpolation points. Here the value of the interpolation function is prescribed, but the value of the Laplacian is not. The corresponding schemes can be derived from the variational form of the applied interpolation. The natural multi-level character of the QT-algorithm makes it simple to define multi-level methods on QT-cell systems. To construct a multi-level method, a sequence of nested cell systems (coarse and fine grids) and proper inter-grid transfer operators (restrictions and prolongations) are required. The cell system of the kth level is defined by simply neglecting the cells the subdivision level of which is greater that k: The restrictions and prolongations can be defined by obvious averaging and weighting techniques, respectively. See Ga´spa´r [11] for details. The essential feature is that though a multi-level QT-based solver is, strictly speaking, a domain type method, the solution procedure is completely determined and controlled by the point set {x1 ; x2 ; …; xN } and does not assume any structure of these points. In this sense, the method is a meshless one.

3. Construction of meshless boundary methods Now let us return to the model problem (1) (Dirichlet problem for the 2D Poisson equation). The effect and the treatment of the boundary and the source term f can be clearly distinguished. First we investigate the effect of the boundary and construct boundary interpolation methods. 3.1. Homogeneous problems The model problem is further simplified to the Laplace equation: Du ¼ 0 in V;

ul›V ¼ u0

ð19Þ

Suppose that the boundary is discretized by the points x1 ; x2 ; …; xN [ ›V: No mesh (grid) structure is assumed. Strictly speaking, Eq. (19) itself can also be considered a harmonic boundary interpolation problem. Indeed,

C. Ga´spa´r / Engineering Analysis with Boundary Elements 28 (2004) 1207–1216

1213

the solution u is harmonic everywhere in V except the interpolation (i.e. the boundary) points and satisfies the interpolation (i.e. the boundary) conditions at these points. However, the number of the ‘interpolation’ points is always infinite. Thus, the singularities arising at the interpolation points (cf. Fig. 1) form a classical single layer potential which exhibits much weaker singularity than the discrete points sources. From numerical point of view, this means that the harmonic interpolation of Section 2.3.1 can be used for solving Dirichlet problems of the Laplace equation only in the case when the number of boundary interpolation points is extremely large (determined by the applied finest gridsize or QT-cellsize, respectively). If, in order to solve the harmonic interpolation problem, a traditional structured grid is used, then we arrive at the classical finite difference method without any computational benefit. In contrast to this, the use of QTs directly to the Laplace (or even to the Poisson) equation is a much more economical technique as pointed out in [8]. However, the boundary must be very densely discretized by the boundary points even if the geometry of the boundary is simple and the solution is smooth. Figs. 3 and 4 show a simple example to illustrate this phenomenon. Here V0 is the unit square, V , V0 is a circle with radius 1/4. The exact solution is uðx; yÞ U x þ y and the boundary condition is consistent with this solution. When a boundary is discretized by 32 points only, logarithmic singularities are generated at the boundary points causing serious error in the approximate solution (the relative error was 13.36% measured in the discrete L2 -norm). Increasing the number of boundary points to 1024, the boundary singularities are cancelled and the relative L2 -error was reduced to 0.10%. The use of the more sophisticated methods based on higher order partial differential operators make it possible to significantly reduce the number of boundary interpolation points. For instance, the biharmonic interpolation based on the boundary points x1 ; x2 ; …; xN produces a smooth function which satisfies the boundary conditions much more exactly than in the case of the harmonic interpolation,

The scaling parameter c should be adjusted in such a way that the ‘computational support’ of the function K0 ðcrÞ appearing in the corresponding fundamental solution is comparable to the characteristic distance of the boundary points. In this case, the interpolation of the boundary condition is still good, while inside the domain V (apart from a narrow neighborhood of the boundary ›V) the interpolation function is nearly harmonic. This implies that the quasi-harmonic boundary interpolation function is a good approximation of the solution of the investigated Dirichlet problem. In our test case, keeping the value

Fig. 3. Harmonic boundary interpolation with boundary singularities.

Fig. 5. Biharmonic boundary interpolation.

Fig. 4. Harmonic boundary interpolation without boundary singularities.

but generally does not satisfy the Laplace equation in V (see Fig. 5), so that the error of approximation is again high (in the presented test case the relative L2 -error was 28.67%; N was set to 32). An optimal compromise can be achieved by using the mixed Laplace– Helmholtz (‘quasi-harmonic’) boundary interpolation. In this approach, the solution of the Dirichlet problem (19) is approximated by the solution of the interpolation problem: DðD 2 c2 IÞu ¼ 0 in V

ð20Þ

uðxk Þ ¼ u0 ðxk Þ ðk ¼ 1; 2; …; NÞ

C. Ga´spa´r / Engineering Analysis with Boundary Elements 28 (2004) 1207–1216

1214

tetraharmonic interpolation: D4 u ¼ 0 in V0 \{x1 ; x2 ; …; xMþN } The boundary conditions along ›V is approximated by the following interpolation conditions: uðxk Þ U u0 ðxk Þ

ðk ¼ M þ 1; M þ 2; …; M þ NÞ;

while the Poisson equation is taken into account also as special interpolation conditions: Duðxk Þ U fk

as described in Section 2.3.3. Note that the variational form of this problem is based on the following subspace of H04 ðV0 Þ :

Fig. 6. Quasi-harmonic boundary interpolation.

N ¼ 32 and setting c to 110, the relative L2 -error was reduced below 0.5% (see Fig. 6). Table 2 shows how depends the relative L2 -error of the approximate solution on the scaling parameter c in the presented test case (with N ¼ 32). 3.2. Inhomogeneous problems

ul›V ¼ u0

ð21Þ

Let us have M points x1 ; x2 ; …; xM [ V scattered inside the domain and N additional points xMþ1 ; xMþ2 ; …; xMþN [ ›V scattered on the boundary. In the method of radial basis functions, two general approaches are known to construct meshless methods: Kansa’s method. The solution u isP approximated by NþM RBFs in the following PNþM form: uðxÞ U j¼1 aj Fj ðx 2 xj Þ: Then DuðxÞ ¼ j¼1 aj DFj ðx 2 xj Þ; and the Poisson equation is required at the points x1 ; x2 ; …; xM ; while the boundary condition is required at the points xMþ1 ; xMþ2 ; …; xMþN : This leads the following system of algebraic equations: NþM X

aj DFj ðxk 2 xj Þ ¼ fk

ðk ¼ 1; 2; …; MÞ

j¼1 NþM X

W U {w [ H04 ðV0 Þ : Dwðx1 Þ ¼ Dwðx2 Þ ¼ · · · ¼ DwðxM Þ ¼ 0; wðxMþ1 Þ ¼ wðxMþ2 Þ ¼ · · · ¼ wðxMþN Þ ¼ 0} Method of particular solutions. In this approach, first the source term f is approximated: f ðxÞ U

Now consider the Poisson equation as a model problem: Du ¼ f in V;

ðk ¼ 1; 2; …; MÞ;

MþN X

aj Fj ðx 2 xj Þ

ð22Þ

j¼1

Let Cj be another RBF defined by DCj ¼ Fj (in practical cases, Cj is known in an analytic form). Having solved the interpolation problem (22), define the function: vðxÞ U PMþN j¼1 aj Cj ðx 2 xj Þ: Then v is (approximately) a particular solution of the model Poisson equation, since DvðxÞ U

MþN X j¼1

aj DCj ðx 2 xj Þ ¼

MþN X

aj Fj ðx 2 xj Þ ¼ f ðxÞ

j¼1

However, v generally does not satisfy the boundary conditions. The solution u is therefore expressed as u ¼ v þ w; where w solves the homogeneous problem with the modified boundary conditions: wl›V ¼ u0 2 vl›V : Using multi-elliptic, e.g. biharmonic interpolation instead of the expression (22), we obtain the following meshless solution algorithm: Step 1. Based on the points x1 ; x2 ; …; xMþN and the values f ðx1 Þ; f ðx2 Þ; …; f ðxMþN Þ; perform a biharmonic interpolation for the function f : DDf ¼ 0 in V0 \{x1 ; x2 ; …; xMþN }

aj Fj ðxk 2 xj Þ ¼ u0 ðxk Þ

j¼1

ðk ¼ M þ 1; M þ 2; …; M þ NÞ In the context of the multi-elliptic, e.g. multi-harmonic interpolation context, this approach leads to the following interpolation strategy. Let u be approximated by a Table 2 The relative L2 -errors with different scaling parameters c 0 50 70 90 110 130 150 200 Relative L2 -error (%) 28.67 5.17 2.95 1.41 0.48 1.09 1.95 3.80

f ðxk Þ is prescribed at x1 ; x2 ; …; xMþN The interpolation is preferably carried out on a QT-cell system generated by the points x1 ; x2 ; …; xMþN : Step 2. The interpolation function obtained in Step 1 is now defined on each cell center belonging to the QT-cell system. Now solve the Poisson equation Dv ¼ f on the same QT-cell system in the larger domain V0 : Along ›V0 ; any regular boundary condition, e.g. Dirichlet boundary condition can be imposed. Step 3. Solve the homogeneous problem Dw ¼ 0 supplied with the modified boundary condition wl›V ¼ u0 2 vl›V : as described in Section 3.1 (by applying a mixed

C. Ga´spa´r / Engineering Analysis with Boundary Elements 28 (2004) 1207–1216

Laplace – Helmholtz boundary interpolation). Then the solution of the inhomogeneous problem is expressed as u ¼ v þ w: Note that the algorithm obviously defines a triharmonictype interpolation for v; which is, however, split into a pair of a biharmonic equation and a Poisson equation. The overall computational cost is typically OððM þ NÞ  logðM þ NÞÞ:

[4]

[5]

[6]

4. Conclusions [7]

A type of meshless methods was constructed. The main idea is the use of the iterated elliptic (biharmonic, biHelmholtz, mixed Laplace – Helmholtz) interpolations. These interpolation methods were proved to be special RBF-type methods based on the fundamental solutions of the applied differential operator. However, the explicit use of these RBFs as well as the solution of large and dense algebraic systems are avoided. In addition to this, neither domain nor boundary integrals are evaluated. The appearing higher order partial differential equations are solved by robust QT-based multi-level methods. The boundary-type version of the presented method based on the mixed Laplace –Helmholtz interpolation was applied to construct a meshless boundary method for the Laplace equation. For inhomogeneous problems, tri- and tetraharmonic-type methods were proposed. They were proved to be a natural generalization of Kansa’s method as well as the method of particular solution in the presented context. It was shown that even the partial differential equations can be considered as special interpolation conditions. In the paper, however, the problems with Neumann and mixed boundary conditions were not investigated since they require a slightly different approach.

[14]

Acknowledgements

[15]

This research was partly sponsored by the Hungarian Scientific Research Fund (OTKA) under the project T34652.

[16]

[8]

[9]

[10]

[11]

[12]

[13]

[17]

References [1] Alves CJS, Chen CS, Saler B. The method of fundamental solutions for solving Poisson problems. In: Brebbia CA, Tadeu A, Popov V, editors. International Series on Advances in Boundary Elements. Proceedings of the 24th International Conference on the Boundary Element Method incorporating Meshless Solutio Seminar, vol. 13. Southatmpon, Boston: WitPress; 2002. p. 67–76. [2] Beatson RK, Light WA. Fast evaluation of radial basis functions: methods for two-dimensional polyharmonic splines. IMA J Numer Anal 1997;17:343 –72. [3] Fam GSA, Rashed YF. A study on the source points locations in the method of fundamental solution. In: Brebbia CA, Tadeu A, Popov V,

[18]

[19]

[20]

[21]

1215

editors. International Series on Advances in Boundary Elements. Proceedings of the 24th International Conference on the Boundary Element Method incorporating Meshless Solutio Seminar, Southatmpon, Boston: WitPress; 2002. p. 297– 312. Fasshauer GE. Solving differential equations with radial basis functions: multilevel methods and smoothing. Advances in computational mathematics, vol. 11. Southampton: WIT Press; 1999. p. 139–59. Floater MS, Iske A. Multistep scattered data interpolation using compactlysupported radial basis functions. J Comput Appl Math 1996;73:65 –78. Franke R. Scattered data interpolation: test of some methods. Math Comput 1982;38(157):181–200. Franke C, Schaback R. Solving partial differential equations by collocation using radial basis functions. Appl Math Comput 1998;93: 73–82. Ga´spa´r C. A multigrid based boundary element method without integral equations. Proceedings of the 13th International Conference on Boundary Element Method, Tulsa, Oklahoma, USA; 1991. Ga´spa´r C, Simbierowicz P. Scattered data interpolation using ¨ , Maksimovic C, unstructured grids. In: Gayer J, Starosolszky O editors. HYDROCOMP’92. Proceedings of the International Conference on Computational Methods and Measurements in Hydraulics and Hydrology, Budapest, Hungary, Budapest: Water Resources Research Centre; 1992. p. 131–8. Ga´spa´r C. Multigrid and multipole techniques in the boundary integral equation method. Boundary elements: implementation and analysis of advanced algorithms. Notes on numerical fluid mechanics, vol. 54. Braunschweig: Vieweg; 1996. p. 102–14. Ga´spa´r C. Flow modelling using quadtrees and multigrid techniques. In: Topping BHV, editor. Proceedings of the Third International Conference on Computational Structures Technology, Budapest, Hungary, Edinburgh: Civil-Comp Press; 1996. p. 31–42. Ga´spa´r G. A multipole expansion technique in solving boundary integral equations. Comput Meth Appl Mech Engng 1998;157: 289–97. Ga´spa´r C. Multigrid technique for biharmonic interpolation with application to dual and multiple reciprocity method. Numer Algorithms 1999;21:165– 83. Ga´spa´r C. Multi-level biharmonic and bi-Helmholtz interpolation with application to the boundary element method. Engng Anal Boundary Elements 2000;24(7–8):559–73. Ga´spa´r C. Fast multi-level meshless methods based on the implicit use of radial basis functions. Lecture notes in computational science and engineering, vol. 26. Berlin: Springer; 2002. p. 143– 60. Golberg MA, Chen CS. The theory of radial basis functions applied to the BEM for inhomogeneous partial differential equations. Boundary Element Commun 1994;5(2):57–61. Golberg MA, Chen CS. A bibliography on radial basis function approximation. Boundary Element Commun 1996;7(4):155–63. Golberg MA, Chen CS, Bowman H. Some recent results and proposals for the use of radial basis functions in the BEM. Engng Anal Boundary Elements 1999;23:285–96. Greaves DM, Borthwick AGL. Hierarchical tree-based finiteelement mesh generation. Int J Numer Meth Engng 1999;45(4): 447 –71. Kiss B, Krebsz A. A sparse matrix representation of the H1=2-norm in finite element spaces. In: Topping BHV, editor. Proceedings of the Third International Conference on Computational Structures Technology, Budapest, Hungary, Edinburgh: Civil-Comp Press; 1996. p. 343–8. Kiss B, Krebsz A. On the Schur complement preconditioners. Comput Struct 1999;73(1-5):537–44.

1216

C. Ga´spa´r / Engineering Analysis with Boundary Elements 28 (2004) 1207–1216

[22] Madych WR, Nelson SA. Multivariate interpolation and conditionally positive definite functions II. Math Comput 1990;54:211– 30. [23] Morinishi K. Gridless type solver—generalized finite difference method. Computational fluid dynamics for the 21st century, notes on numerical fluid mechanics, vol. 78. Berlin: Springer; 2001. p. 43– 58. [24] Partridge PW, Brebbia CA. Computer implementation of the BEM dual reciprocity method for the solution of general field equations. Commun Appl Numer Meth 1990;6:83–92.

[25] Rokhlin V. Rapid solution of integral equations of classical potential theory. J Comput Phys 1985;14:187–207. [26] Wendland H. Piecewise polynomial, positive definite and compactly supported radial basis functions of minimal degree. Adv Comp Math 1995;4:389–96. [27] Yiu KFC, Graeves DM, Cruz S, Saalehi A, Borthwick AGL. Quadtree grid generation—information handling. Boundary fitting and CFD applications. Comput Fluids 1996;25(8):759 –69.