- Email: [email protected]
Particular solutions of 3D Helmholtz-type equations using compactly supported radial basis functions
Engineering Analysis with Boundary Elements 24 (2000) 539–547 www.elsevier.com/locate/enganabound
Particular solutions of 3D Helmholtz-type equations using compactly supported radial basis functions M.A. Golberg a,*, C.S. Chen b, M. Ganesh c a
b
2025 University Circle, Las Vegas, NV 89119, USA Department of Mathematical Sciences, University of Nevada, Las Vegas, NV 89154, USA c School of Mathematics, University of New South Wales, Sydney, NSW 2052, Australia
Abstract In this paper we show how to obtain analytic particular solutions for inhomogeneous Helmholtz-type equations in 3D when the source terms are compactly supported radial basis functions (CS-RBFs) (J. Approx. Theory, 93 (1998) 258). Using these particular solutions we demonstrate the solvability of boundary value problems for the inhomogeneous Helmholtz-type equation by approximating the source term by CS-RBFs and solving the resulting homogeneous equation by the method of fundamental solutions. The proposed technique is a truly mesh-free method and is especially attractive for large-scale industrial problems. A numerical example is given which illustrates the efficiency of the proposed method. 䉷 2000 Elsevier Science Ltd. All rights reserved. Keywords: Helmholtz-type equations; Radial basis functions; Fundamental solution; Particular solution
1. Introduction Following the work of Nardini and Brebbia [1], there has been increasing interest in using the dual reciprocity method (DRM) to solve partial differential equations (PDEs) by boundary methods. Using the DRM or the related method of particular solutions (MPS) enables one to obtain ‘boundary-only’ formulations for inhomogeneous, semi-linear and time-dependent problems by eliminating domain integrals, which typically occur in integral equation formulations for solving boundary value problems. The key to doing this is the ability to analytically calculate particular solutions for various linear PDEs Lu f (without the need to satisfy any boundary conditions). If we can compute approximate particular solutions then the original boundary value problem (with usual subtraction of solutions) can be reduced to solving the homogeneous PDEs Lv 0 with boundary conditions now involving approximate particular solutions. Consequently, boundary integral methods or any scheme which is suited mainly for homogeneous boundary value problems can be used. The process of computing a particular solution is usually done by P approximating the inhomogeneous term f by a series Nj1 aj wj and then solving Lc j wj ; 1 ⱕ j ⱕ N; where {wj } is an appropriate set of linearly independent * Corresponding author. E-mail address: [email protected] (M.A. Golberg).
basis functions. Hence, the choice of {wj } is important and the analysis given in Ref. [2] shows that {wj } needs to provide an accurate approximation to f and should be in the form so that Lc j wj can be solved analytically [2]. In DRM the usual choices of {wj } are the radial basis functions. In recent years, various alternatives to the classical 1 ⫹ r radial basis functions have been considered [3]. Among these are the thin plate and higher order radial splines [4], multiquadrics [5] and Wendland’s compactly supported radial basis functions [6,7]. If L D; the Laplacian, then Lc j wj can be obtained by repeated integration. But for other operators such as Helmholtz-type operators L^ D ^ s2 s 僆 R; different approaches are necessary. In Ref. [4] we obtained analytic particular solutions for L^ when {wj } were polyharmonic radial splines. This was done using the annihilator method [8] and the fact that polyharmonic radial splines are the fundamental solutions for the polyharmonic operators. Because the radial splines are globally supported, approximating f by interpolation leads to full matrices, which become increasingly ill-conditioned as the order of the radial splines increases. As noted in Ref. [6], some of these problems can be alleviated if one uses for {wj } the CS-RBFs. In this case the resulting interpolation matrices are sparse and positivedefinite, which tends to improve the conditioning [6,7]. Sparsity is particularly important, hence CS-RBFs are a natural choice for solving three-dimensional (3D) problems.
0955-7997/00/$ - see front matter 䉷 2000 Elsevier Science Ltd. All rights reserved. PII: S0955-799 7(00)00034-5
540
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
However, the first step in this direction is the ability to calculate analytically particular solutions of the related problem. This is carried out by Chen et al. [6,7] for the Laplacian. In this paper we consider extending the results of Chen et al. [6,7] for the general operator L D ⫺ s2 I in 3D where s 2 is any real or complex number and hence includes both L^ : (We use the notation s 2 for convenience, instead of s to denote any complex number.) These operators are important in the development of algorithms for solving time-dependent PDEs as well as for solving problems in scattering theory. The paper is organized as follows. In Section 2 we give a brief discussion of some boundary value problems which lead naturally to the solution of the inhomogeneous Helmholtz-type equation Lu f and the desirability of using CSRBFs to approximate f. In Section 3 we recall some of the basic approximation properties of CS-RBFs and briefly revisit some results in Refs. [6,7] for Laplace’s equation. In Section 4, we show how to obtain smooth particular solutions for Lc j wj when {wj } are CS-RBFs. Although based on our results these particular solutions can be obtained ‘by hand’, for higher order CS-RBFs, the algebra rapidly becomes formidable and the use of symbolic tools such as Maple or Mathematica, to calculate analytic solutions may be desirable for practical computation. In Section 5 we present numerical examples showing how 3D boundary value problems for Lu f can be solved efficiently by combining the particular solutions obtained in Section 4 with the method of fundamental solutions.
2. Time-dependent and Helmholtz-type problems As motivation for our work, we consider the diffusion boundary value problem Du P; t ut P; t ⫹ f P; t; t ⬎ 0; u P; 0 u0 P; P 僆 V;
1
where V is a bounded domain R d with a C 1 boundary. Such boundary value problems may be solved by using fundamental solutions of the heat operator [9]. But this approach leads to an equivalent equation involving both surface and domain integrals if either f 苷 0 or u0 苷 0 [9]. This approach is in general time-consuming and one loses the benefit of the boundary integral method’s fundamental advantage of reduction in dimension of the space problem. To improve the efficiency of integral equation algorithms much recent research has focused on using time-independent fundamental solutions instead [10–12]. Generally, two approaches have been taken: (i) the elimination of the time dependence by using the Laplace transform [10], and (ii) using finite differences for the time derivative ut [11,12]. For approach (i), we define the Laplace transform in t of
u, say U by U P; l
Z∞
e⫺lt u P; t dt
2
0
and applying this to Eq. (1) gives DU P; l ⫺ lU P; l ⫺u0 P ⫹ F P; l V g P; l;
3
where F P; l is the Laplace transform in t of f. If either u0 苷 0; or f 苷 0; then Eq. (3) is an inhomogeneous Helmholtz-type equation with s2 l; a complex number. Transforming boundary conditions as well [13], one arrives at a boundary value problem which can be solved by using boundary integral techniques. Numerical inversion of the transform gives the solution in time. From the standpoint of modern computers, such algorithms are interesting because they are inherently parallel, as the spatial problems can be solved independent of time. In these formulations domain integrals of the form Z D
G P; Q; lu0 Q dv Q;
Z D
G P; Q; lF P; l dv Q 4
appear where G P; Q; l is the fundamental solution of D ⫺ lI: Since the domain integrals are costly to calculate, it is desirable to avoid them [14]. Hence, much work in this area assumes f u0 0 or f 0 and Du0 0 [10]. However, if we can compute a particular solution of Eq. (3) then we could use stable and efficient algorithms for the corresponding homogeneous problem. To this end, Chen et al. [13] approximated the inhomogeneous term in Eq. (3) by thin plate splines and particular solutions were obtained using the results in Ref. [15]. Muleshkov et al. [4] generalized this using approximations by polyharmonic radial splines. Since radial splines are globally supported, the interpolation matrices which arise in the course of approximation of g P; l are full and become increasingly ill-conditioned as the order of splines and the number of interpolation points increases [16]. These problems become increasingly severe in R3 where hundreds of interpolation points may be needed to obtain good accuracy. One possible approach to alleviating these difficulties is to use CS-RBFs which give rise to sparse interpolation matrices and have improved stability properties because these matrices are positive-definite [6,18]. For Laplace’s equation in R2 ; these properties were successfully exploited in Refs. [6,7] to develop DRM algorithms for Poisson’s equation. In Section 4, these results will be generalized to the Helmholtz-type case. As an alternative to the Laplace transform, one can use finite differencing in ‘t’. This is particularly useful for solving Eq. (1) when f depends on u as well as P; t: For example, defining Un P u P; nt; t ⬎ 0; n 0; 1; 2; …; where t is the time step and approximating ut P; nt ⯝
Un P ⫺ Un⫺1 P t
5
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
gives Vn ; the approximation to Un as the solution to DVn P
Vn P ⫺ Vn⫺1 P ⫹ fn P; t
Then 6
7
Thus Vn satisfies the Helmholtz-type equation with s2 1=t; a real number and f gn : Hence, taking V0 u0 ; boundary value problems for Eq. (1) can be reduced to solving a sequence of Helmholtz-type equations. This method was apparently due to Rothe [11] and has recently been analyzed by Chapko adn Kress [11] and Langdon [12] for equations in R 2. Higher order time-stepping methods may also be developed in this fashion. For example, a Crank–Nicholson scheme gives rise to the sequence Vn P ⯝ u P; nt DVn⫹1 P ⫺
2Vn⫹1 P 2V P ⫺ DVn P ⫺ n t t
8
and Lubich and Schneider considered schemes of the form [9] DVn⫹1 P ⫺ dn⫹1 Vn P
N X
dn⫺k Vk P
9
k0
where {dk } are determined by using A-stable methods for ordinary differential equations. Moreover, as shown in Ref. [21], various time-stepping and/or Laplace transform approaches for non-linear diffusion, wave and convectiondiffusion equations also yield algorithms requiring the solution of inhomogeneous Helmholtz-type equations. Hence, we expect the results in this paper to have wide applicability beyond those already discussed. 3. CS-RBFs For a given f 僆 C V ; to find particular solutions to the linear PDE Lu f
10
we proceed in the following fashion. Assume f to be approximated by a function f N of the form fN P
N X
aj wj P;
P 僆 V
11
j1
and then determine, if possible analytically, a particular solution, say u p;N by solving Lup;N P fN P;
P 僆 V :
12
The function u p;N may be regarded as an ‘approximate’ particular solution to Eq. (10). By linearity, it suffices to solve analytically for unknown c j on V ; where Lc j wj ;
1 ⱕ j ⱕ N:
up;N P
N X
aj c j P:
14
j1
where fn P ⬅ f P; nt: Rearranging Eq. (6) gives V P V P ⫺ n⫺1 ⫹ fn P V gn P: DVn P ⫺ n t t
541
13
If L is a radially and translationally invariant operator, such as the Laplacian D or L D ⫺ s2 I; it is particularly convenient to choose wj as radial basis functions, since solving Eq. (13) can be accomplished by solving the ordinary differential equation [4] Lr c r w r;
15
where L r is the radial part of L. To be more precise, recall that an RBF wj is of the form
wj P w 储P ⫺ Pj 储
16
where w : R⫹ ! R is continuous, Pj 僆 Rd and 储·储 is the Euclidean norm in Rd : Then c j in Eq. (13) is given by
c j P c 储P ⫺ P j 储;
17
where c is a solution of Eq. (15). To obtain efficient algorithms, it is important that w be chosen so that Eq. (15) can be solved analytically and that f N be a good approximation to f. Furthermore, one needs to obtain aj ; j 1; …; N; in Eq. (11) in a computationally efficient manner. The standard way of doing this is to choose N arbitrary distinct points {Pj }N1 in V and then obtain a j ; j 1; …; N by interpolation; i.e. set fN Pk f Pk ;
1 ⱕ k ⱕ N;
18
which requires solving the N linear equations N X
aj wj Pk f Pk ;
1 ⱕ k ⱕ N:
19
j1
If the matrix Aw ⬅ w 储Pj ⫺ Pk 储; 1 ⱕ j; k ⱕ N; is invertible, then the algebraic system (19) has a unique solution for every continuous function f. A sufficient condition for Aw to be invertible is that it be positive-definite, or equivalently, that w be a positive-definite function on R⫹ [17,18]. In general, the commonly used RBFs, such as the conical, polyharmonic splines and multiquadrics are not positive-definite and generally a polynomial term must be added to fN in Eq. (11) in order to guarantee the unique solution of the interpolation problem (19). Moreover, since these functions are globally supported, the matrix Aw is full and can be very ill-conditioned for some RBFs such as multiquadrics [16]. These problems can be particularly severe in R3 : The above problems can be avoided if the RBF w is compactly supported. If w is compactly supported, the matrix Aw is sparse and only a few terms in Eq. (19) need to be evaluated. In his paper, Wu [17] made a significant contribution to the development of CS-RBFs. In 1995, Wendland [18] constructed a new class of positive-definite CS-RBFs. For given smoothness and space dimension, Wendland proved that they are of minimal degree and
542
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
Table 1 are:
Table 1 Wendland’s CS-PD-RBFs in 3D
w w w w
d3
(1 ⫺ (1 ⫺ (1 ⫺ (1 ⫺
r)⫹2 r)⫹4 (4r ⫹ 1) r)⫹6 (35r 2 ⫹ 18r ⫹ 3) r)⫹8 (32r 3 ⫹ 25r 2 ⫹ 8r
8 1=2 > ch ; > > > > < ch3=2 ; 储 f ⫺ fN 储 ∞ ⱕ > > ch5=2 ; > > > : 7=2 ch ;
0
⫹ 1)
C C2 C4 C6
unique up to a constant factor. From a computational point of view, these special features are very attractive. As a result, Chen et al. [6,7] considered the use of Wendland’s CS-RBFs for solving Poission’s equation. Our purpose here is to extend their approach to find particular solutions in R3 to Lu Du ⫺ s2 u f :
20
Before doing this, we review some basic properties of Wendland’s CS-RBFs [18]. Theorem 3.1 ((Wendland)). Up to a constant factor, there exists exactly one positive-definite function w on Rd of the form ( g r; 0 ⱕ r ⱕ 1; w r 21 0; r ⬎ 1; where g r is a univariate polynomial with deg g r d=2 ⫹ 3k ⫹ 1 and satisfying w 僆 C 2k : Here d=2 is the usual integer part, i.e. the unique integer n satisfying n ⱕ d=2 ⬍ n ⫹ 1: The functions g r in the theorem can be explicitly constructed and recursion relations for their coefficients can be found in Ref. [18]. For d 3; the first four CSRBFs and their smoothness are given in Table 1, where 1 ⫺ rl⫹ 1 ⫺ rl for 0 ⱕ r ⬍ 1 and zero otherwise. Defining h sup
min 储P ⫺ P j 储;
22
P僆V 1ⱕjⱕN
it was shown by Wendland in Ref. [18] that if V is compact in R d and satisfies a uniform interior cone condition, then the interpolant fN to f satisfies the uniform error estimate. 储 f ⫺ fN 储∞ ⱕ chk⫹1=2
23
provided that f 僆 H R ; s d=2 ⫹ k ⫹ 1=2; where H s Rd is the Sobolev space of functions satisfying Z 2 2 1 ⫹ 储w储 s 兩 f^ w兩 dw ⬍ ∞ 24 d s
d
R
where f^ is the Fourier transform of f. Hence, for d 3; if f 僆 H k⫹2 R3 and g is of degree 3=2 ⫹ 3k ⫹ 1 3k ⫹ 2; then 储 f ⫺ fN 储∞ ⱕ chk⫹1=2 :
25
So the approximation orders of the interpolants for w ’s in
f 僆 H 2 R3 ; f 僆 H 3 R3 ; f 僆 H 4 R3 ;
26
f 僆 H 5 R3 :
For practical computing purposes, it is important to be able to vary the support of w r: Hence, one can introduce a scaling parameter a ⬎ 0 and define the scaled CS-RBF
wa r w r=a
27
with w as in Eq. (21) or in particular as in Table 1. Using wa r enables one to trade off the sparsity of Aw and the quality of approximation of fN : When a is sufficiently small so that each basis function has only one point Pj in its support, then Aw is diagonal, so has maximum sparsity, but the approximation is poor, since fN is a linear combination of sharp spikes. As a increases, the sparsity of Aw decreases, so Eq. (19) is increasingly costly to solve, but generally the approximation to f improves. 4. Particular solutions with CS-RBFs As indicated in Section 3, with the inhomogeneous term f approximated by the scaled CS-RBF interpolant fN ; in order to calculate particular solutions of Lup;N fN ; it suffices to solve the problem Lr c r w a r; where Lr is the radial part of L. Using Eq. (21) this is equivalent to solving 8 > < g r ; 0 ⱕ r ⱕ a: a 28 Lr c r > : 0; r ⬎ a: In Refs. [6,7], when L D; this was done by direct integration. However, when L D ⫺ s 2 I; this is not convenient and it is preferable to solve Eq. (28) directly as shown below. Of course, our procedure can be used for L D as well, which provides an alternative method to deriving the results in Ref. [6]. Using Eq. (27) and representations of CS-RBFs in Table 1, we have g r=a is of the form g r=a 1 ⫺ r=an p r=a; where p is an appropriately chosen polynomial of degree, say k ⱖ 0; so that w is a C2k CS-RBF. Now using the well-known representation for the radial part of the 3D Laplacian, with L D ⫺ s2 I; Eq. (28) becomes 1 d 2 dc r ⫺ s 2c dr r2 dr 8 n > < 1 ⫺ r p r ; 0 ⱕ r ⱕ a; a a 29 > : 0; r ⬎ a;
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
where for r 0; Eq. (29) is to be considered in the limiting case as r ! 0⫹ : If we let w r ; r ⬎ 0; c r r dj w c j 0 lim⫹ ; j 0; 1; 2; dr j r r!0
Using Eqs. (32) and (34) for 0 ⬍ r ⬍ a;
c 00 r
30
r 2 w 00 r ⫺ 2rw 0 r ⫹ 2w r r3 r2 s2 w r ⫹ rs r ⫺ 2rw 0 r ⫹ 2w r : r3
Since w 0 0; using L’Hospital’s rule, we get
then 1 d 2 dc 1 d 2w r dr r dr 2 r2 dr and hence Eq. (29) is equivalent to 8 n > 2 dr2 : 0;
31
lim c 00 r lim⫹
r!0⫹
×
0 ⱕ r ⱕ a;
r!0
2rs2 w r ⫹ r s r ⫹ r 2 s2 w 0 r ⫹ s r ⫹ r s 0 r ⫺ 2rw 00 r 3r 2
32
r ⬎ a:
lim⫹ r!0
The general solution of Eq. (32) can be shown to be of the form ( ⫺sr A e ⫹ B esr ⫹ q r; 0 ⱕ r ⱕ a; 33 w r r ⬎ a; C e⫺sr ⫹ D esr ; where q r is a particular smooth solution of the first equation in Eq. (32). It is important to note that q r can be determined by the method of undetermined coefficients. However, the task is tedious. Using modern computer symbolic ODE solvers such as Maple or Mathematica, q r can be obtained easily. The four constants in Eq. (33) are to be chosen so that u given by Eq. (30) is twice continuously differentiable at r 0; a and hence on 0; ∞: For the required differentiability of c at 0, the following result shows that it is enough to impose the condition that w 0 0:
s2 w 0 r ⫹ s r ⫹ rs 0 r s 2 w 0 0 ⫹ p 0 : 3 3 36
Finally, using Eqs. (35) and (36) lim⫹
r!0
1 d 2 dc r dr r2 dr
3 lim⫹ r!0
lim⫹ r!0
d2 c 2 dc ⫹ lim⫹ 2 dr r!0 r dr
d2 c lim⫹ s2 c r ⫹ s r: dr 2 r!0
A So by Theorem 4.1 and Eq. (33), the twice continuous differentiability of c at 0 holds if A ⫹ B ⫹ q 0 0:
Theorem 4.1. Let w be a solution of Eq. (32) with w 0 0: Then c defined by Eq. (30) is twice continuously differentiable at 0 with
c 0 w 0 0;
c 0 0 0;
Proof. Since w 0 0; from Eq. (30), clearly c is continuous at 0 with c 0 w 0 0: Further for r ⬎ 0; w 0 r w r rw 0 r ⫺ w r ⫺ 2 : r r r2
34
0
rw r ⫺ w r 0; using L’Hospital’s rule Since lim s 0 0; and Eq. (32), we get r!0⫹
w 00 r lim⫹ sw r ⫹ s r 0: lim⫹ c 0 r lim⫹ 2 r!0 r!0 r!0
A e ⫺sa ⫹ B esa ⫹ q a C e⫺sa ⫹ D esa ; ⫺As e⫺sa ⫹ Bs esa ⫹ q 0 a ⫺Cs e⫺sa ⫹ Ds esa :
Further c r satisfies Eq. (29) as lim r!0 ⫹ :
35
37
Further from Eqs. (30), (32) and (33), it is easy to show that w (and hence u) is twice continuously differentiable at r a ⬎ 0; if (
c 00 0 s 2 w 0 0 ⫹ p 0=3:
c 0 r
543
38
The four unknown constants in Eq. (33) need to satisfy only three equations in Eqs. (37) and (38), one of the variables can be chosen arbitrarily. Since in general there are more interpolation points outside the support (i.e. r ⬎ a; for computational efficiency it is natural to set D 0: Consequently, we get 8 A ⫺B ⫹ q 0; > > > > < e⫺sa q 0 a ⫹ sq a ; B⫺ > 2s > > > : C B e2as ⫺ 1 ⫹ q a eas ⫺ q 0:
39
544
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
As we now know w r explicitly, using Eq. (30) it follows from Theorem 4.1 and Eq. (33) that we have a particular solution of Eq. (29) given by 8 s 2B ⫹ q 0 ⫹ q 0 0; r0 > > > > ⫺sr sr > ⫹ B e ⫹ q r < A e ; 0 ⬍ r ⱕ a; c r r > > > ⫺sr > > : Ce ; r ⬎ a: r
Hence q r; 0 ⱕ r ⱕ a in Eq. (33) is given by q r and
40
where A, B and C are as in Eq. (39). Notice that for r ⬎ a and for Re s ⬎ 0 and large, c r in Eq. (40) can be approximated as
d a2 s2 ⫹ 6 ⫺ 4rs2 a ⫹ 3r 2 s2 q r ⫺ : dr s4 a2 Since
q 0 C e ⫺sr e⫺sr B e2as ⫺ 1 ⫹ q a eas ⫺ q 0 r r 0 1 s a⫺r 1 q a ⫹ sq a e ⫹ ⫺ s a ⫹r 2s r 2s e ⫹ q a eas ⫺ q 0
⬇⫺
e⫺sr r
es a⫺r q 0 a ⫹ sq a 2sr
⫹
4 1 6 2 1 ⫺ 2 r ⫺ 4 2 r ⫹ 2 r2 ⫺ 2 2 r3 as s4 a s s a s a
q a ⫺
2 ; 4 s a
dq s4 a2 ⫹ 6 0 ⫺ 4 2 ; dr sa
dq 6 a ⫺ 4 2 ; dr sa using Eq. (39) we get A⫺
e⫺as 3 ⫹ as ⫹ 4as ; s5 a2
B
e⫺as 3 ⫹ as s5 a2
C
3 eas ⫺ e⫺as 3 ⫹ as ⫺ as 4 ⫹ eas : s5 a2
q a es a⫺r ⫺ q 0 e⫺sr r
⬇ 0 for r q a:
4 ; 4 s a
Using the above values in Eq. (40), we get c: 4.1. Examples Below we consider concrete CS-RBFs given by Eq. (27) with f as in Table 1. The corresponding representations of w r in Eq. (33) and hence c in Eq. (40) are described below. Example 1. tion of
Let w a r 1 ⫺ r=a2⫹ : The general solu-
Example 2. Let w a r 1 ⫺ r=a4⫹ 4 r=a ⫹ 1: In this case, as in the previous example we can show that q r ⫺
480 2880 60 1800 1 ⫺ ⫹ r ⫹ ⫺ s4 a2 s2 s6 a3 s6 a4 s8 a5
240 1440 10 300 3 ⫹ 4 4 ⫹r ⫺ 4 3 ⫺ 6 5 ⫹r sa s2 a2 sa sa 2
d2 w r 2 2 ⫺ s w r 1 ⫺ ; a dr2
0ⱕrⱕa
is given by w r ⫺
⫺4a ⫹ rs2 a2 ⫹ 6r ⫺ 2r2 s2 a ⫹ r3 s2 s4 a2
⫹A e⫺sr ⫹ B esr :
⫺ r4
20 120 ⫹ 4 5 s2 a3 s a
⫹
15 5 4 r ⫺ 2 5 r6 : s2 a4 sa
Hence A, B, C and the solution c can now be easily computed using Eqs. (39) and (40). Example
3.
Let
w a r 1 ⫺ r=a6⫹ 35 r 2 =a ⫹
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
545
ware such as Matematica. As before, the coefficients and the solution c can now be easily computed using Eqs. (39) and (40). 5. Numerical results To substantiate our constructive theoretical results, we consider the following test Helmholtz-type problem (with s 2 400 : D ⫺ s2 Iu f ; Fig. 1. Two connected spheres.
18 r=a ⫹ 3: For this CS-RBF, we have q r ⫺
181440 5040 8467200 504000 5040 ⫺ 6 3 ⫹ ⫹ 8 4 ⫹ 6 2 sa s a s10 a6 sa s8 a 5
725760 3 378 37800 12700800 ⫺ 10 7 ⫹ r ⫺ 2 ⫹ 4 2 ⫹ 6 4 ⫺ s a s sa s10 a7 s a ! 2646000 63000 210 529200 ⫺ 8 5 ⫺ 6 3 ⫺ 4 ⫹ 8 6 sa s a sa sa 90720 2520 362880 252000 ⫺ 4 3 ⫺ 8 7 ⫹ 6 4 sa sa sa s6 a5 ! 2520 4233600 63 6300 441000 ⫹ 4 2 ⫹ ⫹ r3 2 2 ⫹ 4 4 ⫺ 6 5 8 6 a a sa s s s a sa ! 35 2116800 10500 88200 ⫺ 2 ⫺ ⫺ 4 3 ⫹ 6 6 8 7 sa sa sa sa
⫹r 2 ⫺
210 7560 352800 ⫺ 4 5 ⫹ 6 6 s2 a3 s a s a ! 21000 30240 210 ⫹ 4 4 ⫺ 6 7 ⫹ 2 2 sa sa s a
⫹r 4 ⫺
⫹r
5
⫹r 6 ⫹r ⫺
7
! 315 105840 22050 4410 525 ⫺ 6 7 ⫺ 4 5 ⫹ 4 6 ⫺ 2 3 s2 a4 s a sa sa sa ! 700 252 11760 1008 ⫺ 2 5 ⫹ 4 6 ⫺ 4 7 2 4 sa s a s a sa ! ! 105 525 2520 210 18 8 ⫺ 2 5 ⫺ 4 7 ⫹r ⫺ 2 7 s a sa s2 a6 s2 a6 sa
35 9 r sa 2 7
It seems that calculataing q r is quite tedious. However, we remark that q r can be produced and converted to Fortran code automatically by using modern symbolic soft-
u g;
in V;
41
on 2V;
42
where, as shown in Fig. 1, the physical domain V are two connected spheres in R 3 i.e. V { x; y; z 僆 R 3 : H x; y; z ⬍ 1} with H x; y; z min { x ⫺
3 2 4 ; x
⫹
3 2 4 }
⫹ y2 ⫹ z2 :
For numerical test purposes, we chose the inhomogeneous term f in Eq. (41) and the boundary function g in Eq. (42) appropriately so that u x; y; z ex⫹y⫹z =s2 is the exact solution of Eq. (41) ⫺ Eq. (42). To compute approximate solutions of Eq. (41) ⫺ Eq. (42) based our approach, we split the problem into two parts. First, we compute an approximate particular solution, say up;N of Eq. (41) (with N interpolation points in V ) using the constructive formula given in the last section (recall Eqs. (14), (17), (19) and (40)). Then an approximate solution u N of Eq. (41) ⫺ Eq. (42) is uN vN ⫹ up;N ; where vN is a solution of D ⫺ s2 IvN 0;
in V;
v N x; y; z ex⫹y⫹z =s2 ⫺ up;N x; y; z;
43 x; y; z 僆 2V: 44
To evaluate the particular solutions, we chose N 400 quasi-random points as interpolation points in V . The quasirandom points were produced by using subroutine sobseq [20], a quasi-Monte Carlo number generator, to ensure the uniform distribution of the interpolation points. To interpolate the inhomogeneous term in Eq. (41), we chose the CSRBF w a r 1 ⫺ r=a2⫹ (recall Example 1 in the last section). For each choice of cut-off compact support Table 2 Non-zero entries for various compact support a
a
nz
Sparseness (%)
0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8
756 3392 10562 21526 36476 54116 73722 92722 110498
0.40 2.10 6.60 13.5 22.8 33.8 46.0 57.9 69.0
546
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
Table 3 Error estimates for various compact support cut-off parameters
a 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8
储u ⫺ uNM储∞(V ) ⫺2
2.2864 × 10 1.6194 × 10 ⫺2 9.7795 × 10 ⫺3 6.451 × 10 ⫺3 4.5061 × 10 ⫺3 3.8208 × 10 ⫺3 3.3006 × 10 ⫺3 2.9089 × 10 ⫺3 2.6375 × 10 ⫺3
CPU-time (s) 01.40 04.43 07.48 10.63 13.71 16.98 20.53 24.08 27.60
parameter a , these lead to the sparse linear system (19) of size 400 × 400. The sparseness of the matrix (with number of non-zero entries nz) for various values of a are as shown in Table 2. We applied the method of fundamental solutions (MFS) with M 100 quasi-random points on the surface of the physical domain and the same number of interpolation points on the surface of the fictitious domain, a sphere with radius 10, to find an approximate solution vM N of (43) ⫺ (44). For details of the MFS, we refer the reader to Ref. [19] and references therein. So the final computable approximate solution of Eq. (41) ⫺ (42) is uM N vM N ⫹ up;N : We used Matlab on a standard workstation (with multiple users) to compute the numerical solution. The approximate 储u ⫺ uM N 储∞ V error (computed at random 500 points in V for various choices of cut-off parameters a and the corresponding computing time are given in Table 3. The accuracy is particularly good for our high wave number 3D test problem (43) ⫺ (44) with relatively small degrees of freedom. By increasing the number of unknowns or using very efficient solvers for the homogeneous problem (43) ⫺ (44) the accuracy of the solution can be increased. However, the main issue in this paper is to compute approximate particular solutions using CS-RBFs. The test case experiments are in particular encouraging to solve a large class of time-dependent and other applied problems based on our particular solution representation. We intend to do this in our future work.
6. Conclusions We have shown how to compute analytic particular solutions for inhomogeneous Helmholtz-type equations in 3D using CS-RBFs as source terms. It is then shown how to use these particular solutions to develop algorithms for solving boundary value problems for 3D inhomogeneous Helmholtz-type equations by combining these with boundary methods for solving homogeneous Helmholtz-type equations. To facilitate the calculations, symbolic methods are used to obtain the particular solutions. To illustrate the effectiveness of our approach, we solved a boundary value problem for an inhomogeneous Helm-
holtz-type equation by approximating the source term by CS-RBF interpolants and solving the resulting homogeneous equation by the method of fundamental solutions. A study was also made on the effect of the sparsity of the interpolation matrix on CPU time. It was shown that an approximate linear relation holds between CPU time and the support radius of the CS-RBF. As noted, these results are expected to be useful in the development of algorithms for the efficient solution of boundary value problems for a variety of time-dependent partial differential equations. To overcome the difficulties of choosing the cut-off parameter and further improve the efficiency of our approach, the multilevel scheme, a widely used technique in the finite element method, is currently under investigation by our research group. It is expected that the performance of CSRBFs will be further enhanced by this technique.
References [1] Nardini D, Brebbia CA. A new approach to free vibration analysis using boundary elements. Boundary element methods in engineering. Southampton/Berlin/New York: Computational Mechanics/Springer, 1982 (p. 312–26). [2] Golberg MA, Chen CS, Bowman H, Power H. Some comments on the use of radial basis functions in the dual receiprocity method. Comput Mech 1998;21:141–8. [3] Golberg MA, Chen CS, Bowman H. Some recent results and proposals for the use of radial basis functions in the BEM. Engng Anal Bound Elem 1999;23:285–96. [4] Muleshkov AS, Golberg MA, Chen CS. Particular solutions of Helmholtz-type operators using higher order polyharmonic splines. Comput Mech 1999;23:411–9. [5] Golberg MA, Chen CS, Karur SR. Improved multiquadric approximation for partial differential equations. Engng Anal Bound Elem 1996;18:9–17. [6] Chen CS, Brebbia CA, Power H. Dual reciprocity method using compactly supported radial basis functions. Commun Numer Meth Engng 1999;15:137–50. [7] Chen CS, Marcozzi M, Choi S. The method of fundamental solutions and compactly supported radial basis functions — a meshless approach to 3D problems. In: Brebbia CA, Power H, editors. Boundary element methods XXI, Boston/Southampton: WIT Press, 1999. p. 561–70. [8] Golberg MA, Chen CS, Rashed YF. The annihilator method for computing particular solutions to partial differential equations. Engng Anal Bound Elem 1999;23:275–9. [9] Lubich Ch, Schneider R. Time discretization of parabolic boundary integral equations. Numer Math 1992;63:455–81. [10] Zhu S, Satravaha P, Lu X. Solving linear diffusion equations with the dual reciprocity method in Laplace space. Engng Anal Bound Elem 1994;13:1–10. [11] Chapko R, Kress R. Rothe’s method for the heat equation and boundary integral equations. J Int Eq Appl 1997;9:47–68. [12] Langdon S. Domain embedding boundary integral equation methods and parabolic PDEs. PhD thesis, University of Bath, 1999. [13] Chen CS, Rashed YF, Golberg MA. A mesh-free method for linear diffusion equations. Numer Heat Transf Part B 1998;33:469–86. [14] Partridge PW, Brebbia CA, Wrobel LC. The dual reciprocity boundary element method. Southampton/London/New York: Computational Mechanics/Elsevier, 1992. [15] Chen CS, Rashed YF. Evaluation of thin plate spline based particular
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
[16] [17] [18]
[19]
solutions for Helmholtz-type operators for the DRM. Mech Res Commun 1998;25:195–201. Schaback R. Error estimates and condition numbers for radial basis function interpolation. Adv Comput Math 1995;3:251–64. Wu Z. Multivariate compactly supported postive definite radial basis functions. Adv Comput Math 1995;4:283–92. Wendland H. Error estimates for interpolation by compactly supported radial basis functions of minimal degree. J Approx Theory 1998;93:258–73. Golberg MA, Chen CS. The method of fundamental solutions for
547
potential Helmholtz and diffusion problems. Boundary integral methods: numerical and mathematical aspects. Boston/Southampton: WIT Press, 1998. [20] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in Fortran: the art of scientifc computing. Cambridge/New York: Cambridge University Press, 1996. [21] Golberg MA, Chen CS, Muleshkov A. The method of fundamental solutions for time dependent problems. In: Chen CS, Brebbia CA, Pepper DW, editors. Boundary element technology XIII, Boston/ Southampton: WIT Press, 1999. p. 377–86.
Particular solutions of 3D Helmholtz-type equations using compactly supported radial basis functions M.A. Golberg a,*, C.S. Chen b, M. Ganesh c a
b
2025 University Circle, Las Vegas, NV 89119, USA Department of Mathematical Sciences, University of Nevada, Las Vegas, NV 89154, USA c School of Mathematics, University of New South Wales, Sydney, NSW 2052, Australia
Abstract In this paper we show how to obtain analytic particular solutions for inhomogeneous Helmholtz-type equations in 3D when the source terms are compactly supported radial basis functions (CS-RBFs) (J. Approx. Theory, 93 (1998) 258). Using these particular solutions we demonstrate the solvability of boundary value problems for the inhomogeneous Helmholtz-type equation by approximating the source term by CS-RBFs and solving the resulting homogeneous equation by the method of fundamental solutions. The proposed technique is a truly mesh-free method and is especially attractive for large-scale industrial problems. A numerical example is given which illustrates the efficiency of the proposed method. 䉷 2000 Elsevier Science Ltd. All rights reserved. Keywords: Helmholtz-type equations; Radial basis functions; Fundamental solution; Particular solution
1. Introduction Following the work of Nardini and Brebbia [1], there has been increasing interest in using the dual reciprocity method (DRM) to solve partial differential equations (PDEs) by boundary methods. Using the DRM or the related method of particular solutions (MPS) enables one to obtain ‘boundary-only’ formulations for inhomogeneous, semi-linear and time-dependent problems by eliminating domain integrals, which typically occur in integral equation formulations for solving boundary value problems. The key to doing this is the ability to analytically calculate particular solutions for various linear PDEs Lu f (without the need to satisfy any boundary conditions). If we can compute approximate particular solutions then the original boundary value problem (with usual subtraction of solutions) can be reduced to solving the homogeneous PDEs Lv 0 with boundary conditions now involving approximate particular solutions. Consequently, boundary integral methods or any scheme which is suited mainly for homogeneous boundary value problems can be used. The process of computing a particular solution is usually done by P approximating the inhomogeneous term f by a series Nj1 aj wj and then solving Lc j wj ; 1 ⱕ j ⱕ N; where {wj } is an appropriate set of linearly independent * Corresponding author. E-mail address: [email protected] (M.A. Golberg).
basis functions. Hence, the choice of {wj } is important and the analysis given in Ref. [2] shows that {wj } needs to provide an accurate approximation to f and should be in the form so that Lc j wj can be solved analytically [2]. In DRM the usual choices of {wj } are the radial basis functions. In recent years, various alternatives to the classical 1 ⫹ r radial basis functions have been considered [3]. Among these are the thin plate and higher order radial splines [4], multiquadrics [5] and Wendland’s compactly supported radial basis functions [6,7]. If L D; the Laplacian, then Lc j wj can be obtained by repeated integration. But for other operators such as Helmholtz-type operators L^ D ^ s2 s 僆 R; different approaches are necessary. In Ref. [4] we obtained analytic particular solutions for L^ when {wj } were polyharmonic radial splines. This was done using the annihilator method [8] and the fact that polyharmonic radial splines are the fundamental solutions for the polyharmonic operators. Because the radial splines are globally supported, approximating f by interpolation leads to full matrices, which become increasingly ill-conditioned as the order of the radial splines increases. As noted in Ref. [6], some of these problems can be alleviated if one uses for {wj } the CS-RBFs. In this case the resulting interpolation matrices are sparse and positivedefinite, which tends to improve the conditioning [6,7]. Sparsity is particularly important, hence CS-RBFs are a natural choice for solving three-dimensional (3D) problems.
0955-7997/00/$ - see front matter 䉷 2000 Elsevier Science Ltd. All rights reserved. PII: S0955-799 7(00)00034-5
540
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
However, the first step in this direction is the ability to calculate analytically particular solutions of the related problem. This is carried out by Chen et al. [6,7] for the Laplacian. In this paper we consider extending the results of Chen et al. [6,7] for the general operator L D ⫺ s2 I in 3D where s 2 is any real or complex number and hence includes both L^ : (We use the notation s 2 for convenience, instead of s to denote any complex number.) These operators are important in the development of algorithms for solving time-dependent PDEs as well as for solving problems in scattering theory. The paper is organized as follows. In Section 2 we give a brief discussion of some boundary value problems which lead naturally to the solution of the inhomogeneous Helmholtz-type equation Lu f and the desirability of using CSRBFs to approximate f. In Section 3 we recall some of the basic approximation properties of CS-RBFs and briefly revisit some results in Refs. [6,7] for Laplace’s equation. In Section 4, we show how to obtain smooth particular solutions for Lc j wj when {wj } are CS-RBFs. Although based on our results these particular solutions can be obtained ‘by hand’, for higher order CS-RBFs, the algebra rapidly becomes formidable and the use of symbolic tools such as Maple or Mathematica, to calculate analytic solutions may be desirable for practical computation. In Section 5 we present numerical examples showing how 3D boundary value problems for Lu f can be solved efficiently by combining the particular solutions obtained in Section 4 with the method of fundamental solutions.
2. Time-dependent and Helmholtz-type problems As motivation for our work, we consider the diffusion boundary value problem Du P; t ut P; t ⫹ f P; t; t ⬎ 0; u P; 0 u0 P; P 僆 V;
1
where V is a bounded domain R d with a C 1 boundary. Such boundary value problems may be solved by using fundamental solutions of the heat operator [9]. But this approach leads to an equivalent equation involving both surface and domain integrals if either f 苷 0 or u0 苷 0 [9]. This approach is in general time-consuming and one loses the benefit of the boundary integral method’s fundamental advantage of reduction in dimension of the space problem. To improve the efficiency of integral equation algorithms much recent research has focused on using time-independent fundamental solutions instead [10–12]. Generally, two approaches have been taken: (i) the elimination of the time dependence by using the Laplace transform [10], and (ii) using finite differences for the time derivative ut [11,12]. For approach (i), we define the Laplace transform in t of
u, say U by U P; l
Z∞
e⫺lt u P; t dt
2
0
and applying this to Eq. (1) gives DU P; l ⫺ lU P; l ⫺u0 P ⫹ F P; l V g P; l;
3
where F P; l is the Laplace transform in t of f. If either u0 苷 0; or f 苷 0; then Eq. (3) is an inhomogeneous Helmholtz-type equation with s2 l; a complex number. Transforming boundary conditions as well [13], one arrives at a boundary value problem which can be solved by using boundary integral techniques. Numerical inversion of the transform gives the solution in time. From the standpoint of modern computers, such algorithms are interesting because they are inherently parallel, as the spatial problems can be solved independent of time. In these formulations domain integrals of the form Z D
G P; Q; lu0 Q dv Q;
Z D
G P; Q; lF P; l dv Q 4
appear where G P; Q; l is the fundamental solution of D ⫺ lI: Since the domain integrals are costly to calculate, it is desirable to avoid them [14]. Hence, much work in this area assumes f u0 0 or f 0 and Du0 0 [10]. However, if we can compute a particular solution of Eq. (3) then we could use stable and efficient algorithms for the corresponding homogeneous problem. To this end, Chen et al. [13] approximated the inhomogeneous term in Eq. (3) by thin plate splines and particular solutions were obtained using the results in Ref. [15]. Muleshkov et al. [4] generalized this using approximations by polyharmonic radial splines. Since radial splines are globally supported, the interpolation matrices which arise in the course of approximation of g P; l are full and become increasingly ill-conditioned as the order of splines and the number of interpolation points increases [16]. These problems become increasingly severe in R3 where hundreds of interpolation points may be needed to obtain good accuracy. One possible approach to alleviating these difficulties is to use CS-RBFs which give rise to sparse interpolation matrices and have improved stability properties because these matrices are positive-definite [6,18]. For Laplace’s equation in R2 ; these properties were successfully exploited in Refs. [6,7] to develop DRM algorithms for Poisson’s equation. In Section 4, these results will be generalized to the Helmholtz-type case. As an alternative to the Laplace transform, one can use finite differencing in ‘t’. This is particularly useful for solving Eq. (1) when f depends on u as well as P; t: For example, defining Un P u P; nt; t ⬎ 0; n 0; 1; 2; …; where t is the time step and approximating ut P; nt ⯝
Un P ⫺ Un⫺1 P t
5
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
gives Vn ; the approximation to Un as the solution to DVn P
Vn P ⫺ Vn⫺1 P ⫹ fn P; t
Then 6
7
Thus Vn satisfies the Helmholtz-type equation with s2 1=t; a real number and f gn : Hence, taking V0 u0 ; boundary value problems for Eq. (1) can be reduced to solving a sequence of Helmholtz-type equations. This method was apparently due to Rothe [11] and has recently been analyzed by Chapko adn Kress [11] and Langdon [12] for equations in R 2. Higher order time-stepping methods may also be developed in this fashion. For example, a Crank–Nicholson scheme gives rise to the sequence Vn P ⯝ u P; nt DVn⫹1 P ⫺
2Vn⫹1 P 2V P ⫺ DVn P ⫺ n t t
8
and Lubich and Schneider considered schemes of the form [9] DVn⫹1 P ⫺ dn⫹1 Vn P
N X
dn⫺k Vk P
9
k0
where {dk } are determined by using A-stable methods for ordinary differential equations. Moreover, as shown in Ref. [21], various time-stepping and/or Laplace transform approaches for non-linear diffusion, wave and convectiondiffusion equations also yield algorithms requiring the solution of inhomogeneous Helmholtz-type equations. Hence, we expect the results in this paper to have wide applicability beyond those already discussed. 3. CS-RBFs For a given f 僆 C V ; to find particular solutions to the linear PDE Lu f
10
we proceed in the following fashion. Assume f to be approximated by a function f N of the form fN P
N X
aj wj P;
P 僆 V
11
j1
and then determine, if possible analytically, a particular solution, say u p;N by solving Lup;N P fN P;
P 僆 V :
12
The function u p;N may be regarded as an ‘approximate’ particular solution to Eq. (10). By linearity, it suffices to solve analytically for unknown c j on V ; where Lc j wj ;
1 ⱕ j ⱕ N:
up;N P
N X
aj c j P:
14
j1
where fn P ⬅ f P; nt: Rearranging Eq. (6) gives V P V P ⫺ n⫺1 ⫹ fn P V gn P: DVn P ⫺ n t t
541
13
If L is a radially and translationally invariant operator, such as the Laplacian D or L D ⫺ s2 I; it is particularly convenient to choose wj as radial basis functions, since solving Eq. (13) can be accomplished by solving the ordinary differential equation [4] Lr c r w r;
15
where L r is the radial part of L. To be more precise, recall that an RBF wj is of the form
wj P w 储P ⫺ Pj 储
16
where w : R⫹ ! R is continuous, Pj 僆 Rd and 储·储 is the Euclidean norm in Rd : Then c j in Eq. (13) is given by
c j P c 储P ⫺ P j 储;
17
where c is a solution of Eq. (15). To obtain efficient algorithms, it is important that w be chosen so that Eq. (15) can be solved analytically and that f N be a good approximation to f. Furthermore, one needs to obtain aj ; j 1; …; N; in Eq. (11) in a computationally efficient manner. The standard way of doing this is to choose N arbitrary distinct points {Pj }N1 in V and then obtain a j ; j 1; …; N by interpolation; i.e. set fN Pk f Pk ;
1 ⱕ k ⱕ N;
18
which requires solving the N linear equations N X
aj wj Pk f Pk ;
1 ⱕ k ⱕ N:
19
j1
If the matrix Aw ⬅ w 储Pj ⫺ Pk 储; 1 ⱕ j; k ⱕ N; is invertible, then the algebraic system (19) has a unique solution for every continuous function f. A sufficient condition for Aw to be invertible is that it be positive-definite, or equivalently, that w be a positive-definite function on R⫹ [17,18]. In general, the commonly used RBFs, such as the conical, polyharmonic splines and multiquadrics are not positive-definite and generally a polynomial term must be added to fN in Eq. (11) in order to guarantee the unique solution of the interpolation problem (19). Moreover, since these functions are globally supported, the matrix Aw is full and can be very ill-conditioned for some RBFs such as multiquadrics [16]. These problems can be particularly severe in R3 : The above problems can be avoided if the RBF w is compactly supported. If w is compactly supported, the matrix Aw is sparse and only a few terms in Eq. (19) need to be evaluated. In his paper, Wu [17] made a significant contribution to the development of CS-RBFs. In 1995, Wendland [18] constructed a new class of positive-definite CS-RBFs. For given smoothness and space dimension, Wendland proved that they are of minimal degree and
542
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
Table 1 are:
Table 1 Wendland’s CS-PD-RBFs in 3D
w w w w
d3
(1 ⫺ (1 ⫺ (1 ⫺ (1 ⫺
r)⫹2 r)⫹4 (4r ⫹ 1) r)⫹6 (35r 2 ⫹ 18r ⫹ 3) r)⫹8 (32r 3 ⫹ 25r 2 ⫹ 8r
8 1=2 > ch ; > > > > < ch3=2 ; 储 f ⫺ fN 储 ∞ ⱕ > > ch5=2 ; > > > : 7=2 ch ;
0
⫹ 1)
C C2 C4 C6
unique up to a constant factor. From a computational point of view, these special features are very attractive. As a result, Chen et al. [6,7] considered the use of Wendland’s CS-RBFs for solving Poission’s equation. Our purpose here is to extend their approach to find particular solutions in R3 to Lu Du ⫺ s2 u f :
20
Before doing this, we review some basic properties of Wendland’s CS-RBFs [18]. Theorem 3.1 ((Wendland)). Up to a constant factor, there exists exactly one positive-definite function w on Rd of the form ( g r; 0 ⱕ r ⱕ 1; w r 21 0; r ⬎ 1; where g r is a univariate polynomial with deg g r d=2 ⫹ 3k ⫹ 1 and satisfying w 僆 C 2k : Here d=2 is the usual integer part, i.e. the unique integer n satisfying n ⱕ d=2 ⬍ n ⫹ 1: The functions g r in the theorem can be explicitly constructed and recursion relations for their coefficients can be found in Ref. [18]. For d 3; the first four CSRBFs and their smoothness are given in Table 1, where 1 ⫺ rl⫹ 1 ⫺ rl for 0 ⱕ r ⬍ 1 and zero otherwise. Defining h sup
min 储P ⫺ P j 储;
22
P僆V 1ⱕjⱕN
it was shown by Wendland in Ref. [18] that if V is compact in R d and satisfies a uniform interior cone condition, then the interpolant fN to f satisfies the uniform error estimate. 储 f ⫺ fN 储∞ ⱕ chk⫹1=2
23
provided that f 僆 H R ; s d=2 ⫹ k ⫹ 1=2; where H s Rd is the Sobolev space of functions satisfying Z 2 2 1 ⫹ 储w储 s 兩 f^ w兩 dw ⬍ ∞ 24 d s
d
R
where f^ is the Fourier transform of f. Hence, for d 3; if f 僆 H k⫹2 R3 and g is of degree 3=2 ⫹ 3k ⫹ 1 3k ⫹ 2; then 储 f ⫺ fN 储∞ ⱕ chk⫹1=2 :
25
So the approximation orders of the interpolants for w ’s in
f 僆 H 2 R3 ; f 僆 H 3 R3 ; f 僆 H 4 R3 ;
26
f 僆 H 5 R3 :
For practical computing purposes, it is important to be able to vary the support of w r: Hence, one can introduce a scaling parameter a ⬎ 0 and define the scaled CS-RBF
wa r w r=a
27
with w as in Eq. (21) or in particular as in Table 1. Using wa r enables one to trade off the sparsity of Aw and the quality of approximation of fN : When a is sufficiently small so that each basis function has only one point Pj in its support, then Aw is diagonal, so has maximum sparsity, but the approximation is poor, since fN is a linear combination of sharp spikes. As a increases, the sparsity of Aw decreases, so Eq. (19) is increasingly costly to solve, but generally the approximation to f improves. 4. Particular solutions with CS-RBFs As indicated in Section 3, with the inhomogeneous term f approximated by the scaled CS-RBF interpolant fN ; in order to calculate particular solutions of Lup;N fN ; it suffices to solve the problem Lr c r w a r; where Lr is the radial part of L. Using Eq. (21) this is equivalent to solving 8 > < g r ; 0 ⱕ r ⱕ a: a 28 Lr c r > : 0; r ⬎ a: In Refs. [6,7], when L D; this was done by direct integration. However, when L D ⫺ s 2 I; this is not convenient and it is preferable to solve Eq. (28) directly as shown below. Of course, our procedure can be used for L D as well, which provides an alternative method to deriving the results in Ref. [6]. Using Eq. (27) and representations of CS-RBFs in Table 1, we have g r=a is of the form g r=a 1 ⫺ r=an p r=a; where p is an appropriately chosen polynomial of degree, say k ⱖ 0; so that w is a C2k CS-RBF. Now using the well-known representation for the radial part of the 3D Laplacian, with L D ⫺ s2 I; Eq. (28) becomes 1 d 2 dc r ⫺ s 2c dr r2 dr 8 n > < 1 ⫺ r p r ; 0 ⱕ r ⱕ a; a a 29 > : 0; r ⬎ a;
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
where for r 0; Eq. (29) is to be considered in the limiting case as r ! 0⫹ : If we let w r ; r ⬎ 0; c r r dj w c j 0 lim⫹ ; j 0; 1; 2; dr j r r!0
Using Eqs. (32) and (34) for 0 ⬍ r ⬍ a;
c 00 r
30
r 2 w 00 r ⫺ 2rw 0 r ⫹ 2w r r3 r2 s2 w r ⫹ rs r ⫺ 2rw 0 r ⫹ 2w r : r3
Since w 0 0; using L’Hospital’s rule, we get
then 1 d 2 dc 1 d 2w r dr r dr 2 r2 dr and hence Eq. (29) is equivalent to 8 n > 2
31
lim c 00 r lim⫹
r!0⫹
×
0 ⱕ r ⱕ a;
r!0
2rs2 w r ⫹ r s r ⫹ r 2 s2 w 0 r ⫹ s r ⫹ r s 0 r ⫺ 2rw 00 r 3r 2
32
r ⬎ a:
lim⫹ r!0
The general solution of Eq. (32) can be shown to be of the form ( ⫺sr A e ⫹ B esr ⫹ q r; 0 ⱕ r ⱕ a; 33 w r r ⬎ a; C e⫺sr ⫹ D esr ; where q r is a particular smooth solution of the first equation in Eq. (32). It is important to note that q r can be determined by the method of undetermined coefficients. However, the task is tedious. Using modern computer symbolic ODE solvers such as Maple or Mathematica, q r can be obtained easily. The four constants in Eq. (33) are to be chosen so that u given by Eq. (30) is twice continuously differentiable at r 0; a and hence on 0; ∞: For the required differentiability of c at 0, the following result shows that it is enough to impose the condition that w 0 0:
s2 w 0 r ⫹ s r ⫹ rs 0 r s 2 w 0 0 ⫹ p 0 : 3 3 36
Finally, using Eqs. (35) and (36) lim⫹
r!0
1 d 2 dc r dr r2 dr
3 lim⫹ r!0
lim⫹ r!0
d2 c 2 dc ⫹ lim⫹ 2 dr r!0 r dr
d2 c lim⫹ s2 c r ⫹ s r: dr 2 r!0
A So by Theorem 4.1 and Eq. (33), the twice continuous differentiability of c at 0 holds if A ⫹ B ⫹ q 0 0:
Theorem 4.1. Let w be a solution of Eq. (32) with w 0 0: Then c defined by Eq. (30) is twice continuously differentiable at 0 with
c 0 w 0 0;
c 0 0 0;
Proof. Since w 0 0; from Eq. (30), clearly c is continuous at 0 with c 0 w 0 0: Further for r ⬎ 0; w 0 r w r rw 0 r ⫺ w r ⫺ 2 : r r r2
34
0
rw r ⫺ w r 0; using L’Hospital’s rule Since lim s 0 0; and Eq. (32), we get r!0⫹
w 00 r lim⫹ sw r ⫹ s r 0: lim⫹ c 0 r lim⫹ 2 r!0 r!0 r!0
A e ⫺sa ⫹ B esa ⫹ q a C e⫺sa ⫹ D esa ; ⫺As e⫺sa ⫹ Bs esa ⫹ q 0 a ⫺Cs e⫺sa ⫹ Ds esa :
Further c r satisfies Eq. (29) as lim r!0 ⫹ :
35
37
Further from Eqs. (30), (32) and (33), it is easy to show that w (and hence u) is twice continuously differentiable at r a ⬎ 0; if (
c 00 0 s 2 w 0 0 ⫹ p 0=3:
c 0 r
543
38
The four unknown constants in Eq. (33) need to satisfy only three equations in Eqs. (37) and (38), one of the variables can be chosen arbitrarily. Since in general there are more interpolation points outside the support (i.e. r ⬎ a; for computational efficiency it is natural to set D 0: Consequently, we get 8 A ⫺B ⫹ q 0; > > > > < e⫺sa q 0 a ⫹ sq a ; B⫺ > 2s > > > : C B e2as ⫺ 1 ⫹ q a eas ⫺ q 0:
39
544
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
As we now know w r explicitly, using Eq. (30) it follows from Theorem 4.1 and Eq. (33) that we have a particular solution of Eq. (29) given by 8 s 2B ⫹ q 0 ⫹ q 0 0; r0 > > > > ⫺sr sr > ⫹ B e ⫹ q r < A e ; 0 ⬍ r ⱕ a; c r r > > > ⫺sr > > : Ce ; r ⬎ a: r
Hence q r; 0 ⱕ r ⱕ a in Eq. (33) is given by q r and
40
where A, B and C are as in Eq. (39). Notice that for r ⬎ a and for Re s ⬎ 0 and large, c r in Eq. (40) can be approximated as
d a2 s2 ⫹ 6 ⫺ 4rs2 a ⫹ 3r 2 s2 q r ⫺ : dr s4 a2 Since
q 0 C e ⫺sr e⫺sr B e2as ⫺ 1 ⫹ q a eas ⫺ q 0 r r 0 1 s a⫺r 1 q a ⫹ sq a e ⫹ ⫺ s a ⫹r 2s r 2s e ⫹ q a eas ⫺ q 0
⬇⫺
e⫺sr r
es a⫺r q 0 a ⫹ sq a 2sr
⫹
4 1 6 2 1 ⫺ 2 r ⫺ 4 2 r ⫹ 2 r2 ⫺ 2 2 r3 as s4 a s s a s a
q a ⫺
2 ; 4 s a
dq s4 a2 ⫹ 6 0 ⫺ 4 2 ; dr sa
dq 6 a ⫺ 4 2 ; dr sa using Eq. (39) we get A⫺
e⫺as 3 ⫹ as ⫹ 4as ; s5 a2
B
e⫺as 3 ⫹ as s5 a2
C
3 eas ⫺ e⫺as 3 ⫹ as ⫺ as 4 ⫹ eas : s5 a2
q a es a⫺r ⫺ q 0 e⫺sr r
⬇ 0 for r q a:
4 ; 4 s a
Using the above values in Eq. (40), we get c: 4.1. Examples Below we consider concrete CS-RBFs given by Eq. (27) with f as in Table 1. The corresponding representations of w r in Eq. (33) and hence c in Eq. (40) are described below. Example 1. tion of
Let w a r 1 ⫺ r=a2⫹ : The general solu-
Example 2. Let w a r 1 ⫺ r=a4⫹ 4 r=a ⫹ 1: In this case, as in the previous example we can show that q r ⫺
480 2880 60 1800 1 ⫺ ⫹ r ⫹ ⫺ s4 a2 s2 s6 a3 s6 a4 s8 a5
240 1440 10 300 3 ⫹ 4 4 ⫹r ⫺ 4 3 ⫺ 6 5 ⫹r sa s2 a2 sa sa 2
d2 w r 2 2 ⫺ s w r 1 ⫺ ; a dr2
0ⱕrⱕa
is given by w r ⫺
⫺4a ⫹ rs2 a2 ⫹ 6r ⫺ 2r2 s2 a ⫹ r3 s2 s4 a2
⫹A e⫺sr ⫹ B esr :
⫺ r4
20 120 ⫹ 4 5 s2 a3 s a
⫹
15 5 4 r ⫺ 2 5 r6 : s2 a4 sa
Hence A, B, C and the solution c can now be easily computed using Eqs. (39) and (40). Example
3.
Let
w a r 1 ⫺ r=a6⫹ 35 r 2 =a ⫹
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
545
ware such as Matematica. As before, the coefficients and the solution c can now be easily computed using Eqs. (39) and (40). 5. Numerical results To substantiate our constructive theoretical results, we consider the following test Helmholtz-type problem (with s 2 400 : D ⫺ s2 Iu f ; Fig. 1. Two connected spheres.
18 r=a ⫹ 3: For this CS-RBF, we have q r ⫺
181440 5040 8467200 504000 5040 ⫺ 6 3 ⫹ ⫹ 8 4 ⫹ 6 2 sa s a s10 a6 sa s8 a 5
725760 3 378 37800 12700800 ⫺ 10 7 ⫹ r ⫺ 2 ⫹ 4 2 ⫹ 6 4 ⫺ s a s sa s10 a7 s a ! 2646000 63000 210 529200 ⫺ 8 5 ⫺ 6 3 ⫺ 4 ⫹ 8 6 sa s a sa sa 90720 2520 362880 252000 ⫺ 4 3 ⫺ 8 7 ⫹ 6 4 sa sa sa s6 a5 ! 2520 4233600 63 6300 441000 ⫹ 4 2 ⫹ ⫹ r3 2 2 ⫹ 4 4 ⫺ 6 5 8 6 a a sa s s s a sa ! 35 2116800 10500 88200 ⫺ 2 ⫺ ⫺ 4 3 ⫹ 6 6 8 7 sa sa sa sa
⫹r 2 ⫺
210 7560 352800 ⫺ 4 5 ⫹ 6 6 s2 a3 s a s a ! 21000 30240 210 ⫹ 4 4 ⫺ 6 7 ⫹ 2 2 sa sa s a
⫹r 4 ⫺
⫹r
5
⫹r 6 ⫹r ⫺
7
! 315 105840 22050 4410 525 ⫺ 6 7 ⫺ 4 5 ⫹ 4 6 ⫺ 2 3 s2 a4 s a sa sa sa ! 700 252 11760 1008 ⫺ 2 5 ⫹ 4 6 ⫺ 4 7 2 4 sa s a s a sa ! ! 105 525 2520 210 18 8 ⫺ 2 5 ⫺ 4 7 ⫹r ⫺ 2 7 s a sa s2 a6 s2 a6 sa
35 9 r sa 2 7
It seems that calculataing q r is quite tedious. However, we remark that q r can be produced and converted to Fortran code automatically by using modern symbolic soft-
u g;
in V;
41
on 2V;
42
where, as shown in Fig. 1, the physical domain V are two connected spheres in R 3 i.e. V { x; y; z 僆 R 3 : H x; y; z ⬍ 1} with H x; y; z min { x ⫺
3 2 4 ; x
⫹
3 2 4 }
⫹ y2 ⫹ z2 :
For numerical test purposes, we chose the inhomogeneous term f in Eq. (41) and the boundary function g in Eq. (42) appropriately so that u x; y; z ex⫹y⫹z =s2 is the exact solution of Eq. (41) ⫺ Eq. (42). To compute approximate solutions of Eq. (41) ⫺ Eq. (42) based our approach, we split the problem into two parts. First, we compute an approximate particular solution, say up;N of Eq. (41) (with N interpolation points in V ) using the constructive formula given in the last section (recall Eqs. (14), (17), (19) and (40)). Then an approximate solution u N of Eq. (41) ⫺ Eq. (42) is uN vN ⫹ up;N ; where vN is a solution of D ⫺ s2 IvN 0;
in V;
v N x; y; z ex⫹y⫹z =s2 ⫺ up;N x; y; z;
43 x; y; z 僆 2V: 44
To evaluate the particular solutions, we chose N 400 quasi-random points as interpolation points in V . The quasirandom points were produced by using subroutine sobseq [20], a quasi-Monte Carlo number generator, to ensure the uniform distribution of the interpolation points. To interpolate the inhomogeneous term in Eq. (41), we chose the CSRBF w a r 1 ⫺ r=a2⫹ (recall Example 1 in the last section). For each choice of cut-off compact support Table 2 Non-zero entries for various compact support a
a
nz
Sparseness (%)
0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8
756 3392 10562 21526 36476 54116 73722 92722 110498
0.40 2.10 6.60 13.5 22.8 33.8 46.0 57.9 69.0
546
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
Table 3 Error estimates for various compact support cut-off parameters
a 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8
储u ⫺ uNM储∞(V ) ⫺2
2.2864 × 10 1.6194 × 10 ⫺2 9.7795 × 10 ⫺3 6.451 × 10 ⫺3 4.5061 × 10 ⫺3 3.8208 × 10 ⫺3 3.3006 × 10 ⫺3 2.9089 × 10 ⫺3 2.6375 × 10 ⫺3
CPU-time (s) 01.40 04.43 07.48 10.63 13.71 16.98 20.53 24.08 27.60
parameter a , these lead to the sparse linear system (19) of size 400 × 400. The sparseness of the matrix (with number of non-zero entries nz) for various values of a are as shown in Table 2. We applied the method of fundamental solutions (MFS) with M 100 quasi-random points on the surface of the physical domain and the same number of interpolation points on the surface of the fictitious domain, a sphere with radius 10, to find an approximate solution vM N of (43) ⫺ (44). For details of the MFS, we refer the reader to Ref. [19] and references therein. So the final computable approximate solution of Eq. (41) ⫺ (42) is uM N vM N ⫹ up;N : We used Matlab on a standard workstation (with multiple users) to compute the numerical solution. The approximate 储u ⫺ uM N 储∞ V error (computed at random 500 points in V for various choices of cut-off parameters a and the corresponding computing time are given in Table 3. The accuracy is particularly good for our high wave number 3D test problem (43) ⫺ (44) with relatively small degrees of freedom. By increasing the number of unknowns or using very efficient solvers for the homogeneous problem (43) ⫺ (44) the accuracy of the solution can be increased. However, the main issue in this paper is to compute approximate particular solutions using CS-RBFs. The test case experiments are in particular encouraging to solve a large class of time-dependent and other applied problems based on our particular solution representation. We intend to do this in our future work.
6. Conclusions We have shown how to compute analytic particular solutions for inhomogeneous Helmholtz-type equations in 3D using CS-RBFs as source terms. It is then shown how to use these particular solutions to develop algorithms for solving boundary value problems for 3D inhomogeneous Helmholtz-type equations by combining these with boundary methods for solving homogeneous Helmholtz-type equations. To facilitate the calculations, symbolic methods are used to obtain the particular solutions. To illustrate the effectiveness of our approach, we solved a boundary value problem for an inhomogeneous Helm-
holtz-type equation by approximating the source term by CS-RBF interpolants and solving the resulting homogeneous equation by the method of fundamental solutions. A study was also made on the effect of the sparsity of the interpolation matrix on CPU time. It was shown that an approximate linear relation holds between CPU time and the support radius of the CS-RBF. As noted, these results are expected to be useful in the development of algorithms for the efficient solution of boundary value problems for a variety of time-dependent partial differential equations. To overcome the difficulties of choosing the cut-off parameter and further improve the efficiency of our approach, the multilevel scheme, a widely used technique in the finite element method, is currently under investigation by our research group. It is expected that the performance of CSRBFs will be further enhanced by this technique.
References [1] Nardini D, Brebbia CA. A new approach to free vibration analysis using boundary elements. Boundary element methods in engineering. Southampton/Berlin/New York: Computational Mechanics/Springer, 1982 (p. 312–26). [2] Golberg MA, Chen CS, Bowman H, Power H. Some comments on the use of radial basis functions in the dual receiprocity method. Comput Mech 1998;21:141–8. [3] Golberg MA, Chen CS, Bowman H. Some recent results and proposals for the use of radial basis functions in the BEM. Engng Anal Bound Elem 1999;23:285–96. [4] Muleshkov AS, Golberg MA, Chen CS. Particular solutions of Helmholtz-type operators using higher order polyharmonic splines. Comput Mech 1999;23:411–9. [5] Golberg MA, Chen CS, Karur SR. Improved multiquadric approximation for partial differential equations. Engng Anal Bound Elem 1996;18:9–17. [6] Chen CS, Brebbia CA, Power H. Dual reciprocity method using compactly supported radial basis functions. Commun Numer Meth Engng 1999;15:137–50. [7] Chen CS, Marcozzi M, Choi S. The method of fundamental solutions and compactly supported radial basis functions — a meshless approach to 3D problems. In: Brebbia CA, Power H, editors. Boundary element methods XXI, Boston/Southampton: WIT Press, 1999. p. 561–70. [8] Golberg MA, Chen CS, Rashed YF. The annihilator method for computing particular solutions to partial differential equations. Engng Anal Bound Elem 1999;23:275–9. [9] Lubich Ch, Schneider R. Time discretization of parabolic boundary integral equations. Numer Math 1992;63:455–81. [10] Zhu S, Satravaha P, Lu X. Solving linear diffusion equations with the dual reciprocity method in Laplace space. Engng Anal Bound Elem 1994;13:1–10. [11] Chapko R, Kress R. Rothe’s method for the heat equation and boundary integral equations. J Int Eq Appl 1997;9:47–68. [12] Langdon S. Domain embedding boundary integral equation methods and parabolic PDEs. PhD thesis, University of Bath, 1999. [13] Chen CS, Rashed YF, Golberg MA. A mesh-free method for linear diffusion equations. Numer Heat Transf Part B 1998;33:469–86. [14] Partridge PW, Brebbia CA, Wrobel LC. The dual reciprocity boundary element method. Southampton/London/New York: Computational Mechanics/Elsevier, 1992. [15] Chen CS, Rashed YF. Evaluation of thin plate spline based particular
M.A. Golberg et al. / Engineering Analysis with Boundary Elements 24 (2000) 539–547
[16] [17] [18]
[19]
solutions for Helmholtz-type operators for the DRM. Mech Res Commun 1998;25:195–201. Schaback R. Error estimates and condition numbers for radial basis function interpolation. Adv Comput Math 1995;3:251–64. Wu Z. Multivariate compactly supported postive definite radial basis functions. Adv Comput Math 1995;4:283–92. Wendland H. Error estimates for interpolation by compactly supported radial basis functions of minimal degree. J Approx Theory 1998;93:258–73. Golberg MA, Chen CS. The method of fundamental solutions for
547
potential Helmholtz and diffusion problems. Boundary integral methods: numerical and mathematical aspects. Boston/Southampton: WIT Press, 1998. [20] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in Fortran: the art of scientifc computing. Cambridge/New York: Cambridge University Press, 1996. [21] Golberg MA, Chen CS, Muleshkov A. The method of fundamental solutions for time dependent problems. In: Chen CS, Brebbia CA, Pepper DW, editors. Boundary element technology XIII, Boston/ Southampton: WIT Press, 1999. p. 377–86.