Recent developments in the numerical evaluation of particular solutions in the boundary element method

NOIfl'H. H ( K I A N D

Recent D e v e l o p m e n t s in the Numerical Evaluation of Particular Solutions in the Boundary Element M e t h o d Michael Golberg

2025 University Circle Las Vegas, Nevada 89119-6051

Transmitted by John Casti

ABSTRACT We survey recent work on the numerical evaluation of particular solutions in the boundary elements method (BEM), particularly rbf approximations and related mathematical convergence theory. Current research suggests, at least in R 2, that the best choice for bases is the generalized multiquadrics. Numerical results are presented to support this observation.

1.

INTRODUCTION

In recent years there has been considerable interest in computing particular solutions for inhomogeneous partial differential equations. T h a t research, most often associated with the Dual and Multiple Reciprocity Methods (DRM and MRM), has resulted in numerous papers [1] and at least two books [2, 3]. Despite that, there are still a considerable number of unresolved theoretical and computational issues, particularly those related to convergence and stability [4]. In this respect a number of results have recently appeared, which go some way toward resolving these problems [5-7].

APPLIED MATHEMATICS AND COMPUTATION 75:91-101 (1996) © Elsevier Science Inc., 1996 0096-3003/96/$15.00 655 Avenue of the Americas, New York, NY 10010 SSDI 0096-3003(95)00123-Y

92

M. GOLBERG

Much of the work in this area through 1990 is discussed in the book by Partridge, Brebbia, and Wrobel [2], so here we concentrate primarily on work since then and present some new results we and others have obtained. Although a number of approaches have been taken, here we focus on the use of radial basis function (rbf) approximations, currently the preferred method

[2]. In Section 2 we review some of the theory of: i'bfs as it~applies :to :the approximation of particular solutions for rotationally symmetric differential operators. In particular, we discuss some recent convergenceresults for rbf interpolants given by Schmidt and Strese [6], Madych [5], Schaback and Wu [4] and Powell [7]. These theoretical results along with a number of experimental studies [5, 8] suggest that multiquadrics axe the best choice for basis functions. Numerical results given in Section 4 bear this out. In Section 3 we examine the convergence of the DRM using rbfs and we conclude in Section 5 with some directions for future research. 2.

RADIAL BASIS FUNCTION APPROXIMATIONS

2.1. Computational Considerations We consider boundary value problems (bvps) for equations of the form

Lu( P) = f ( P ) ,

P ~ D,

(2.1)

where D is a domain in R d, d = 2, 3 with boundary S and L is an elliptic differential operator. A particular solution to (2.1) is a solution up which does not necessarily satisfy the given boundary conditions. To solve bvps for (2.1) it is convenient to reduce it to an equivalent homogeneous equation, as the resulting boundary integral reformulation will be free of domain integrals [2]. Hence, let u = v + up where Lv = 0 with boundary conditions for v derived from u and up. The function v can be determined by a standard boundary elements method (BEM); the remaining problem is to determine up. If L has constant coefficients and f is a polynomial or trigonometric function, then up can be determined by the method undetermined coefficients [2, 9]. More generally, one can represent up as a volume potential

up(P) = foG( P, Q)f( Q) dV( Q),

(2.2)

Numerical Evaluation of Particular Solutions

93

where G(P, Q) is a fundamental solution for L. Usually (2.2) must be evaluated numerically, which may require complicated meshing of D and this is what we wanted to avoid in the first place. In [9] Atkinson showed that replacing D with D, an ellipse containing D in R 2 or an ellipsoid in R 3, and L = A, the Laplacian, that (2.2) could be evaluated using standard numerical integration methods without explicit discretization of D. This approach can be used for any constant coefficient operator. In [10, 11] we showed, using rbf approximations to f, that the domain integrals in Atkinson's formula could be reduced to line or surface integrals as in the DRM. Another possibility, which we have yet to implement, is to reduce Atkinson's formula to a series of boundary integrals as in the MRM. To avoid the direct computation of domain integrals Nardini and Brebbia in [12] proposed the following method, now known as the DRM: assume that f can be approximated sufficiently accurately by a finite series of basis functions (hi, i = 1 , 2 , . . . , N; i.e., N

f ( P ) ~ ] ( P ) = E A~¢i(P),

(2.3)

i=1

then replace f by ] in (1) and approximate up by ~ where N

L~ -

E A~bi-

(2.4)

i=1

By linearity, ~t = F~N=1h3b, where L~b, = ~b~,

i = 1 , 2 , . . . , N.

(2.5)

In [12] and in almost all subsequent papers on the DRM [1, 2], it was required that (2.5) be solved analytically. If L had constant coefficients, choosing {~bi} to be polynomials or trigonometric functions is reasonable, as {~b~} may be obtained by undetermined coefficients. However, finding such approximations to f can be tricky, particularly if collocation is used to obtain f. In contrast to the one-dimensional case, due to a theorem of Haar [13], one generally cannot interpolate with polynomial or trigonometric bases, as unisolvent sets of functions must depend on the interpolation points themselves. This problem occurred in the paper of Cheng et al. [14], where the singularity of the interpolation matrix was observed without noting the cause.

94

M. GOLBERG

On the other hand, polynomial approximations can be obtained by Taylor expansions as in [9], and trigonometric approximations by Fourier series. The latter can be computationally costly [9, 15], but can be alleviated to some extent by the fast Fourier transform as shown by Kassab and Nordlund in [15]. Because interpolation is usually the most efficient method of obtaining approximations to f, it is important to have basis functions for which the interpolation problem is uniquely solvable. An important way of doing this is to choose {4)i} to be rbfs; 4)i(P) = 4 ) ( J I P - QilI), where 4): R+-* R is continuous, {Qi }N i= a c_ D U S is a set of interpolation points and IJ" JJ is the Euclidean norm on R a. Although this choice does not appear to be intuitively obvious, in their original paper on the DRM Nardini and Brebbia chose 4)(r) = 1 + r [12]. Until recently, virtually all papers on the DRM have used these bases without further comment. As we observed in [11, 16], one way of motivating rbfs as "natural" bases is via a theorem by Duchon [17], who showed that the thin plate splines (TPS) 4)(r) - r 2 log r in R 2 and 4)(r) = r in R 3 solved the optimal interpolation problem--hence is the natural generalization of cubic splines. As is well known, the function f1 which interpolates to f at { xi}iN 1 c R and minimizes f~_=]f"(x)J 2 dx is a cubic spline. In R d, d = 2, 3, the analog of this problem is to determine fl which interpolates to f at {Q,}~ 1 ({ Qi}ig l are noncollinear in R 2 and nonplanar in R a) and minimizes

I( I)

fJ,R a( i , E/ ) = l[OxiO~xjJ dx.

(2.6)

QJl)f,(P) = F-N=1hir~ log r i = EN=laiyi + a x + b y + c, E ,N= l h i ~__ E i =Nl h i X i _ = 0 and in R a ( P = ( z , y, z))h(P) = 2ff=iA, ri + ax + by + cy + d, 2ff=lhi = Y'.~=lhiXi = ]~ff=~hz y~ = Z~= ~hi z~ = 0. Thus the optimal interpolant is the sum of an rbf As shown in [17] in R ~ ( P = ( x, y), r i -- IIP -

and a polynomial. An important byproduct of Duchon's theorem is the unique solvability of the interpolation problem for TPS, which follows from the uniqueness of the optimal interpolant. In [11], using TPS along with the method of fundamental solutions (MFS), we solved a number of bvps with improved accuracy over the DRM with 1 + r bases. Despite their optimal properties, TPS have a number of drawbacks. First they are only C 1 in R 2 and are nondifferentiable in R a, a property they share with the 1 + r rbfs. This has proved to be a problem in a number of applications of the DRM [2, 18]. Moreover, recent analysis has shown a rather slow rate of convergence. For instance, Powell has shown for f

Numerical Evaluation of Particular Solutions

95

C2(D), and ft the T P S interpolant on D c R 2 that IIf - fzlL • chllog hi where h = SUppe omin{ll P - Qill} [7]. In R a Madych and Nelson [19] and Wu and Schaback [4] have given local error bounds which are 0(hl/2). Similar results hold for the 1 + r functions. One possible remedy for these problems is to use higher order Duchon splines [17]. These are interpolants which minimize the functional

amf In(y) = frt ~'~ ax~ . . . . ~x~ a

]2

(2.7)

dx

where the sum in (2.7) is taken over all multi-indices (O~1, . . . , O/d) satisfying la~ m. As shown by Duchon in [17] the optimal interpolant in this case is of the form f I ( P ) = ~=lh,tbm(ri) +Pm where d~m(r) = r 2~-d, d odd, and thin(r) = r 2~- d log r, d even, and Pm is a polynomial degree m - 1 satisfying ( t = d i m P ~ , the set of polynomials of degree ~< m - 1) ZN=IAjbi(Qj) = 0, i = 1 , . . . , t, where {bi}~=1 is a basis for Pro" Unfortunately, the implementation of these higher order splines is awkward as the collocation points must he chosen as part of a unisolvent set for Pro. However, Duchon's theorem indicates that rbf approximations a r e " natural" rather than ad hoc as some authors have suggested [2, 20]. Even though splines are optimal in the previous sense, current research suggests that the best choice for rbf interpolants for smoothness and accuracy is the family of multiquadrics d ~ ( r ) = ( r 2 + c~) ~/2 where fl is an odd integer and c is a parameter to be determined by the data. For fl = 1 / 2 these interpolants were introduced by Hardy in 1971 to interpolate geophysical data [21] and have been intensely studied since then. An excellent review of this literature through 1990 is given in [22]. In 1982, Franke compared 29 two-dimensional interpolation methods available at the time and found that the Hardy multiquadrics performed best, followed by T P S [8]. However, it was not until Michelli's landmark 1986 paper that their mathematical properties began to be understood [23]. In that paper it was shown that ~bl/2(r) was conditionally positive definite of order zero showing that the interpolation problem =

N

z0 +

~ Ai~bl/2(llQi i=1

N

-

Qjll) = f(Qj),

j = 1 , 2 , . . . , N, ~ A, = 0, i=l

(2.8)

96

M. GOLBERG

was uniquely solvable for arbitrary sets of interpolation points { Qj}N 1 in general position. Since this result, the convergence and stability properties of the multiquadrics have been actively investigated [5, 24]. In particular, Madych and Nelson [18] and Wu and Schaback [4] have recently shown that multiquadric interpolants converge exponentially fast in h, and Madych has shown that for certain classes of functions (the properties are given in terms of the Fourier transform of f ) that convergence is exponential in c [5]. Similar properties hold for the generalized multiquadrics ~bg(r). Motivated by these results, we have begun to investigate the application of generalized multiquadric bases in the BEM and MFS. In particular, we have shown by appropriate choice of c, that one can achieve several orders of magnitude improvement in accuracy over the 1 + r and TPS bases. Further improvement can be obtained by proper placement of the collocation points. An outstanding problem in this area is the "optimal" choice of ( fl, c). If L is the two-dimensional Laplacian, we have found that /3 = - 3 is a good choice because in this case the basis functions are positive definite [5] and the particular solutions computed from (2.5) are analytically simpler than for/3 = 1. The inverse multiquadrics ~b l(r) also work well and have been used by Ingber and Phan-Tien [25] and Ingber [26] with c = 1 in the DRM for the heat equation. The generalized multiquadrics for/3 = - 3 , - 1, 1 and c = 1 have also been proposed by Zheng et al. in [27]--but one should note that the particular solutions given there are in error. On the other hand, choosing c seems to be more difficult. This problem has been investigated in the geophysics literature, and for that type of data a number of heuristics and ad-hoc formulas have been proposed. For some f's, as mentioned previously, the error decreases exponentially in c, so large c's should be optimal. However, as c increases the bases become "flatter" and the conditioning of the interpolation matrix deteriorates, thus increasing round-off error [5]. We have found this type of behavior in numerical experiments. In [28] Kansa showed that further improvement in accuracy can be obtained by choosing c to depend on Q~. Some theoretical justification for this is given in [29].

2.2. Computing Particular Solutions Traditionally in the DRM once ] has been obtained ~ is calculated by (2.4) and (2.5). For rbf approximations f( P) = E~=lA~¢(lIP- Q~rl) + p where p is a polynomial. By linearity, ~ = ~i=lAi~bi N + X where Lq*i = ¢i and LX = p. If L has constant coefficients, then X can be obtained by the method of undetermined coefficients. In general, the analytic solution to

Numerical Evaluation of Particular Solutions

97

(2.5) is feasible only if L is rotationally symmetric. In practice, the two most common operators of this type are the Laplacian and the biharmonic operator A2. If L is rotationally symmetric, then it is natural to look for radial solutions O(r). In this case L~b(r)= L(r, d / dr)~b( r), where L is an ordinary differential operator depending only on r. Thus (2.5) becomes L~b(r) = ~b(r), which can be integrated by standard variation of parameters formulas. Care should be taken to add an appropriate solution to L~b = 0 so that $ is continuous at r = 0. Specific examples for A, A2 and the Helmholtz operators A + k2 can be found in [2, 30]. When L is not rotationally symmetric, other approaches need to be taken to calculate ~. For example, if L = a232/3x 2 + b2~2/~ y2, then choosing "radial" approximations to match the symmetry of L may be appropriate. Letting r A = ( x 2 / a 2 + y2/b2)l/2, L ~ ( r ~ ) = g(rA) , SO expanding f in functions depending on r A would seem natural. For general operators where one cannot integrate (2.5) analytically, one can look for direct approximations to up as follows. Assume ~ = y,N=1high, + p and determine {h i} and p by solving L~(Qj) -= O, i = 1, 2 , . . . , N along with the consistency equations ~Nj=IAjb,(Qj)= 0, i = 1 , 2 , . . . , t. In [20] Allesandri and Tralli used this approach to determine bicubic spline approximations ~, and Kansa [28] has used this technique to solve the full bvp as an alternative to the finite element or finite difference methods. This method is under investigation.

3.

CONVERGENCE OF T H E DRM W I T H rbfs

Despite recent advances in understanding the convergence of rbf approximations, this does not immediately imply convergence of ~, even if obtained analytically from ~ To see the difficulty, recall that particular solutions are not unique. Thus, what particular solution, if any, does ~ converge to? Of course convergence of ~ is not of primary importance; rather convergence of the corresponding algorithm for solving the given bvp is. In [6] Schmidt and Strese apparently gave a proof of the DRMBEM using 1 + r rbfs for Poisson's equation. However, we believe there is a gap in the proof that needs to be filled. We examine their argument to see where the difficulty lies. We consider the bvp (2.1) with L = A and Dirichlet boundary conditions ul s = g. To solve this bvp let f be an rbf approximation to f with satisfying A~---)~ Let v = u - ~ and let vh be an approximation to v obtained by some version of the BEM and u h an approximation to u given

98

M. GOLBERG

by u h = vh + ~ . T h e e r r o r e h = u - u h = v - u h + u - w ( w = v + ~ ) a n d we consider convergence of u h in L2(/9). By the triangle inequality

(3.1)

Ilenll2 < IIv - vhll2 -4-Ilu - wl12.

T o estimate the error IIu - wlf2 we use the well-known, a-priori estimate for solutions to Poisson's equation [31].

fD ( U-- W) 2 dV ~ Clfs( u - w) 2 US -[- c2fD[h ( u -- w)] 2 dV. Now u becomes

w l s = O and A ( u -

w)=Au-Av-AS=f-/so

f n ( u - w) 2 dV <<. c2f~ ( f -

f)2 dV

(3.2)

t h a t (3.2)

(3.3)

and taking square roots in (3.3) gives

Ilu- wll < cll/-/11 .

(3.4)

fli2

F r o m our previous c o m m e n t s ] i f = O(h"), a > 0, for TPS, multiquadric and 1 + r rbfs. Hence, I l u - wl]2 ~ 0 , h ~ 0 in these cases. It remains to bound [[v - vhll2 to prove convergence of u h. In [6] it is claimed t h a t this follows immediately from standard error estimates for the BEM. This a p p a r e n t l y is not the case, because the b o u n d a r y d a t a g - uls depend on ~ which in turn depends on the discretization p a r a m e t e r h. More precisely, suppose the b v p is solved by a projection method with defining projections {Ph}. Then s t a n d a r d error estimates show t h a t [30] IIv - vhll2 ~< cl(h) + ]1~ - Ph~ll where I1" II is an appropriate surface norm on S and Q(h) ~ O, h --* O. For typical discretizations ]]~ - Phil] ~< chl3ll~lls, /3 > 0, where ]]-Ils is a suitable Sobolev norm. In general, we have no information on fibrils. Until such bounds are available, the convergence of the D R M appears to be an open problem. 4.

NUMERICAL RESULTS

Here we give a brief s u m m a r y of some recent numerical results for Poisson's equation using the MFS and various rbf approximations to f. For comparison with previous work we have solved the b v p A u = - - x 2, ul s = 0

Numerical Evaluation of Particular Solutions

99

on the ellipse { x ~ + 4 y~ ~< 1} in R 2. In [2] this problem was solved using the BEM with linear elements and 1 + r basis elements and in [20] with linear elements and bicubic spline interpolants. In [11] we used the MFS with T P S interpolants. Here we consider the use of the MFS and generalized multiquadrics with /3 = - 3 , - 1 , 1 and various values of c. The results in Tables 1 and 2, the solution at (0, 0), show a substantial increase in accuracy (2 orders of magnitude) over previous BEM and MFS results. Our MFS calculations were done with 16 sources evenly spaced around a circle of radius eight to solve the bvp. The approximate particular solutions were obtained using thin plate spline and multiquadric interpolants on the 33 points given in [2]. The exact solution at (0, 0) is 0.1365836. In Table 1 numbers are rounded to three decimal places. In Table 2 the numbers in parentheses are errors. The starred value is the most accurate result. 5.

CONCLUSIONS

We have surveyed recent work on the numerical evaluation of particular solutions in the BEM, particularly rbf approximations, and related mathematical convergence theory. Current research suggests, at least in R 2, that the best choice for bases is the generalized multiquadrics. Numerical results are presented to support this observation. Further 2D and 3D calculations are needed to further substantiate these conclusions.

TABLE 1 MFS SOLUTION OF A u =

-- X2 -- COMPARISON

WITH THE BEM

MFS Solution

DRM Solution

H-bicubic Solution

0.135

0.127

0.132

MFS SOLUTIONOF A u ~-

TABLE 2 X2 WITHMULTIQUADRICS(~(r)

c=l /3 = - 3 /3 = - 1 /3 = 1

0.13745537 (0.87 × 10 -a) 0.136660043 (0.75 × 10 -4) 0.136224666 (-0.36 X 10 -a)

-

-

c=5 0.136583442 (0.19 × 10-5) 0.136581780 (0.36 × 10-5) 0.136581376 (-0.39 × 10 -5)

C=10 0.136588786 (0.34 × 10 -5) 0.136601790 (0.16 × 10 -4) 0.136585078* (0.28 × 10 -6)

100

M. GOLBERG

The author acknowledges the generous help of Professor C. S. Chen in the preparation of this paper. REFERENCES 1 M. A. Golberg, The Numerical Evaluation of Particular Solutions in the BEM--A Review, Boundary Elements Communications, 6:99-106 (1995). 2 P. W. Partridge, L. C. Wrobel, and C. A. Brebbia, The Dual Reciprocity Boundary Element Method, Computational Mechanics Publications, 1991. 3 A. J. Nowak and A. C. Neves, The Multiple Reciprocity Boundary Element Method, Computational Mechanics Publications, 1994. 4 Z. Wu and R. Shaback, Local error estimates for radial basis function interpolation of scattered data, IMA J. Num. Anal. 13:13-27 (1993). 5 W . R . Madych. Miscellaneous error bounds for multiquadric and related interpt> lators, Comp. Math. Appl. 24:121-138 (1992). 6 G. Schmidt and H. Strese. BEM for Poisson equation. Eng. Anal. Boundary Elements 10:119-123 (1992). 7 M. J. D. Powell, The Uniform Convergence of Thin Plate @lines in Two Dimensions, University of Cambridge Numerical Analysis Report DAMTP 1993/NA 16, 1993. 8 R. Franke. Scattered data interpolation: tests for some methods. Math. Comp. 38:181-199 (1982). 9 K. Atkinson. The numerical evaluation of particular solutions for Poisson's equation. IMA J. Num. Anal. 5:319-338 (1985). 10 M.A. Golberg and C. S. Chen. On a method of Atkinson for evaluating domain integrals in the BEM, Appl. Math. Comp. 60:125-138 (1994). 11 M . A . Golberg. The method of fundamental solutions for Poisson's equation, Eng. Anal. Boundary Elements (in press). 12 D. Nardini and C. A. Brebbia, A New Approach to Free Vibration Analysis Using Boundary Elements, Boundary Element Methods in Engineering, Springer-Verlag, 1983. 13 K. Salkanskas, Some Relationshzps Between Surface Splines and Kriging, ISNM61, Birkh~user Verlag Basel, 1982. 14 A . H . Cheng, O. Lafe, and S. Grilli, Dual reciprocity BEM method based on complete set global shape functions, in Boundary Elements X V (C. A. Brebbia and J. Rencis, Eds.), Computational Mechanics Publications, Vol. 1, 1993, pp. 343-357. 15 A . J . Kassab and R. S. Nordlund. Efficient implementation of the Fourier dual reciprocity BEM using two-dimensional fast Fourier transforms, Eng. Anal. Boundary Elements 12:93-102 (1993). 16 M . A . Golberg and C. S. Chen. The theory of radial basis functions applied to the BEM for inhomogeneous partial differential equations, Boundary Elements Communications 5:57-61 (1994). 17 J. Duchon, Splines Minimizing Rotation Invariant Seminorm in Sobolev Spaces, Constructive Theory of Functions of Several Variables, Lecture notes in Mathematics 571 (W. Schempp and K. Zeller, Eds.), Springer-Verlag, 1978.

Numerical Evaluation of Particular Solutions

101

18 L. Lee and T. W. Wu, Applications of the dual reciprocity method to acoustic radiation in a subeanic nonuniform flow, in Boundary Element Technology IX (C. A. Brebbia and A. J. Kassab, eds.), Computational Mechanics Publications, 1994. 19 W . R . Madych and S. A. Nelson. Multivariate interpolation and conditionally positive definite functions II, Math. Comp. 54:211-230 (1990). 20 C . A . Allesandri and A. Tralli. A spline-based approach for avoiding domain integrations in the BEM, Computers and Structures 41:859-868 (1991). 21 R.L. Hardy. Multiquadric equations of topography and other irregular surfaces, J. Geophys. Res. 26:1905-1915 (1971). 22 R. L. Hardy. Theory and applications of the multiquadric-biharmonic method--20 years of discovery 1968-1988, Computers Math. Applic. 19:163-208 (1990). 23 C.A. Micchelli. Interpolation of scattered data: distance matrices and conditionally positive definite functions, Constr. Approx. 2:11-22 (1986). 24 M.D. Buhmann and N. Dyn, Error estimates for multiquadric interpolation, in Curves and Surfaces (P. J. Laurent, A. Le Mehaute, and L. S. Schumaker, Eds.), Academic Press, 1991, pp. 51-58. 25 M.S. Ingber and N. Phan-Tien. A boundary element approach for parabolic differential equations using a class of particular solutions, Appl. Math. Modelling 16:124-132 (1992). 26 M.S. Ingber, A triple reciprocity BEM for transient heat conduction analysis, in Boundary Element Technology IX (C. A. Brebbia and A. J. Kassab, Eds.), Computational Mechanics Publications, 10, 1994. 27 R. Zheng, C. J. Coleman, and N. Phan-Tien. A boundary element approach for nonhomogeneous potential problems, Computational Mechanics 7:279-288 (1991). 28 E. J. Kansa. Multiquadrics--A scattered data approximation scheme with applications to computational fluid dynamics-I, Surface approximations and partial derivative estimates, Computers Math. Applic. 19:127-145 (1990). 29 M. D. Buhmann and C. A. Micchelli. Multiquadric interpolation improved, Computers Math. AppIic., 24:21-26 (1992). 30 S. Zhu, Particular solutions associated with the Helmholtz operator used in the BEM, Boundary Element Abstracts 4:231-233 (1993). 31 K . E . Custafson, Introduction to Partial Differential Equations and Hilbert Space Methods, John Wiley and Sons, New York, 1980.