Proceeding of China Congress on Computational Mechanics’ 2001

RBF-based meshless boundary knot method and boundary
particle method
Wen CHEN*
Department of Informatics, University of Oslo, P.O.Box 1080, Blindern, 0316 Oslo, Norway
(E-mail: wenc@ifi.uio.no)
Classification: meshfree methods
Abstract: This paper is concerned with the two new boundary-type radial basis function collocation schemes,
boundary knot method (BKM) and boundary particle method (BPM). The BKM is developed based on the dual
reciprocity theorem, while the BKM employs the multiple reciprocity technique. Unlike the method of fundamental
solution, the two methods use the non-singular general solution instead of singular fundamental solution to
circumvent the controversial artificial boundary outside physical domain. Compared with the boundary element
method, both BKM and BPM are meshless, super-convergent, meshfree, integration-free, symmetric, and
mathematically simple collocation techniques for general PDE’s. In particular, the BPM does not require any inner
nodes for inhomogeneous problems. In this study, the accuracy and efficiency of the two methods are numerically
demonstrated to some 2D, 3D Helmholtz and convection-diffusion problems under complicated geometries.
Keyword: boundary knot method; boundary particle method; radial basis function; meshfree; method of fundamental
solution, dual reciprocity BEM, multiple reciprocity BEM.

1. Introduction
In last decade much effort has been devoted to developing a variety of meshless schemes for numerical
partial differential equation (PDE) discretization. The driving force behind the scene is that the meshbased methods such as the standard FEM and BEM often require prohibitively computational effort to
mesh or remesh in handling high-dimensional, moving boundary, and complex-shaped boundary
problems. Many of the meshfree techniques available now are based on using moving least square
strategy (MLS). In most cases, a shadow element is still necessary for numerical integration rather than
for function interpolation. Therefore, these methods are not truly meshfree.
Exceptionally, the methods based on radial basis function (RBF) are inherently meshfree due to the
fact that the RBF method does not employ the MLS at all and uses the one-dimensional distance variable
irrespective of dimensionality of problems. Therefore, the RBF methods are independent of
dimensionality and complexity of problem geometry. Nardini and Brebbia[1] in 1982 have actually applied
the RBF concept to develop currently popular dual reciprocity BEM (DR-BEM) without a notion of
“RBF” and the use of then related advances in multivariate scattered data processing. Only after Kansa’s
pioneer work [2] in 1990, the research on the RBF method for PDE’s has become very active.
Among the existing RBF schemes, the so-called Kansa’s method is a domain-type collocation
technique, while the method of fundamental solution (MFS) (regular BEM versus singular BEM) is a
typical boundary-type RBF collocation methodology. The MFS outperforms the standard BEM in terms
of integration free, convergence speed, easy-to-use, and meshfree merits[3]. The main drawback of the
MFS is due to use the fictitious boundary outside physical domain. The arbitrariness in the determination
of the artificial boundary introduces such troublesome issues as stability and accuracy in dealing with
complicated geometry systems and impedes the efficacy of the MFS to practical engineering problems.
Instead of using the singular fundamental solution, Chen and Tanaka[4,5] exploited the non-singular
general solution to the approximation of homogeneous solution and removed the controversial artificial
*

This research is supported by Norwegian Research Council. The final submission is on 27 Aug. 2001.

Proceeding of China Congress on Computational Mechanics’ 2001

boundary in the MFS. The method is called the boundary knot method (BKM). Like the MFS and
DRBEM, the BKM also employs the dual reciprocity method to approximate particular solution. For
recent development of the BKM see ref. [6]. Some preliminary numerical experiments[4-7] show that the
BKM can produce excellent results with relatively a small number of nodes for various linear and
nonlinear problems.
On the other hand, in recent years the multiple reciprocity BEM (MR-BEM)[8] has attracted much
attention in the BKM community. The method uses high-order fundamental solutions to approximate
high-order homogeneous solutions and then get the particular solutions. The advantages and
disadvantages of the MRM relative to the DR-BEM are that it does not use inner nodes at all for
inhomogeneous problems but that it requires more computing effort. By analogy with the MR-BEM,
Chen[6,7] introduced the meshfree boundary particle method (BPM) which combines the RBF and multiple
reciprocity principle to formulate a simple and efficient boundary-only meshless collocation method.
The purpose of this study is to introduce the BKM and BPM and testify them to some typical 2D and
3D PDE systems under complicated geometry. The paper is concluded with some remarks.

2. Boundary knot method
To clearly illustrate our idea, consider the following example without loss of generality
ℜ{u} = f (x ) ,
x∈Ω ,
x ⊂ Su ,
u (x ) = R (x ) ,
∂u ( x )
= N (x ) ,
∂n

x ⊂ ST ,

where ℜ is differential operator, x means multi-dimensional independent variable, and n is the unit
outward normal. The solution of Eq. (1) can be expressed as
u = uh + u p ,
where uh and up are the homogeneous and particular solutions, respectively. The latter satisfies

{ }

ℜ u p = f (x )

(1)
(2a)
(2b)

(3)
(4)

but does not necessarily satisfy boundary conditions.
To evaluate the particular solution, the inhomogeneous term is approximated first by
f (x ) ≅

∑ λ ϕ (r ) ,

N +L

j

j

(5)

j =1

where λ j are the unknown coefficients. N and L are respectively the numbers of knots on the domain and
boundary. The use of interior points is usually necessary to guarantee the accuracy and convergence of
the BKM solution. r j = x − x j represents the Euclidean distance norm, and ϕ is the radial basis function.
By forcing approximation representation (5) to exactly satisfy Eq. (4) at all nodes, we can uniquely
determine
λ = Aϕ−1{ f (xi )} ,
(6)
where Aϕ is nonsingular RBF interpolation matrix. Finally, we can get particular solutions at any point by
summing localized approximate particular solutions
up =

∑λ φ( x − x ),

N +L

j

j

(7)

j =1

where the RBF φ is related to the RBF ϕ through operator ℜ. Substituting Eq. (6) into Eq. (7) yields
u p = ΦAϕ−1 { f (x i )} ,

(8)

where Φ is a known matrix comprised of φ (rij).
On the other hand, the homogeneous solution uh has to satisfy both governing equation and boundary
conditions. Unlike the dual reciprocity BEM[1] and MFS[3] using the singular fundamental solution, the
BKM[4-8] approximates uh by means of nonsingular general solution, namely,

Proceeding of China Congress on Computational Mechanics’ 2001
L

uh (x ) = ∑α k u # (rk ) ,

(9)

k =1

where k is the index of source points on boundary; u# is the nonsingular general solution of operator ℜ. αk
are the desired coefficients. Collocating Eqs. (1) and (2a,b) at all boundary and interior knots in terms of
representation (8) and (9), we have the unsymmetric BKM schemes. For the sake of brevity, the
respective details are omitted here[4-7]. In order to get symmetric BKM scheme for self-adjoint operators,
we modify the BKM approximate expression (9) to homogeneous solution uh as
Ld + L N

Ld

uh (x ) = ∑ asu # (rs ) −

∑

s =1

as

s = L d +1

∂u # (rs )
,
∂n

(10)

where n is the unit outward normal as in boundary condition (2b), and Ld and LN are respectively the
numbers of knots at the Dirichlet and Neumann boundary surfaces. The minus sign associated with the
second term is due to the fact that the Neumann condition of the first order derivative is not self-adjoint.
In terms of expression (10), the collocation analogue equations (1a) and (2a,b) are written as
Ld

Ld + L N

∑ a u (r ) − ∑ a
∂u (r )
∑ a ∂n − ∑
#

s

is

s

s =1

s = L d +1

#

Ld

∂u # (ris )
= R (xi ) − u p (xi ) ,
∂n

( ) = N (x ) − ∂u (x ) ,

∂ 2u # rjs

Ld + L N

js

as

s

s =1

∂n 2

s = L d +1

Ld

∑ asu # (rls ) −
s =1

Ld + LN

∑

as

s = L d +1

(11)
p

j

∂u # (rls )
= ul − u p (xl ) .
∂n

∂n

j

(12)
(13)

The system matrix of the above equations is symmetric if operator ℜ{} is self-adjoint. Note that i, s and j
are reciprocal indices of Dirichlet ( Su ) and Neumann boundary ( SΓ ) nodes. l indicates response knots
inside domain Ω. After the solution of the above simultaneous algebraic equations, we can employ the
obtained expansion coefficients α and inner knot solutions ul to calculate the BKM solution at any knot. It
is stressed here that the MFS could not produce the symmetric interpolation matrix in any way.
The present form of the BKM uses the expansion coefficient rather than the direct physical variable in
the approximation of boundary value. Therefore, such BKM is called as the indirect BKM. Chen et al. [6,7]
also gave the direct BKM with physical variable as basic variable.

3. Boundary particle methods
Just like the DR-BEM, the interior nodes are usually necessary in the BKM for inhomogeneous
problems. A rival to the DR-BEM is the MR-BEM[8], which applies the multiple reciprocity principle to
circumvent the domain integral without using any inner node. In this section, we will develop a boundaryonly RBF scheme called the boundary particle method [6,7] with the multiple reciprocity principle.
The multiple reciprocity method assumes that the particular solution in Eq. (3) can be approximated by
higher-order homogeneous solution, namely,
∞

u = uh0 + u 0p = uh0 + ∑ uhm ,

(14)

m =1

where superscript m is the order index of homogeneous solution. uh0 and u 0p are equivalent to
homogeneous solution uh and particular solution up in Eq. (3). Through an incremental differentiation
operation via operator ℜ{}, we have successively higher order differential equations:
 uh0 (xi ) = R(xi ) − u 0p (xi )
 0
∂u 0p x j ,
 ∂u h x j
N
x
=
−
j

∂n
 ∂n

( ) ( )
 ℜ {u (x )} = f (x ) − ℜ {u (x )}

 ∂ℜ {u (x )} ∂ ( f (x ) − ℜ {u (x )}) ,
=
( )

0

0




1
h

1
h

∂n

0

i

1
p

i

0

j

j

∂n

i

1
p

j

(15a)

(15b)

Proceeding of China Congress on Computational Mechanics’ 2001

{ }
{ ( )} (

{

}
{ ( )}) ,

 ℜ n −1 uhn (xi ) = ℜ n − 2 { f (xi )} − ℜ n −1 u np (xi )
 n −1 n
∂ ℜ n − 2 f x j − ℜ n −1 u np x j
 ∂ℜ u h x j
=

∂n
∂n


{ ( )}

n=2,3,…,

(15c)

where ℜn{} denotes the n-th order operator ℜ{}, say ℜ1{}=ℜℜ0{}, and ℜ0{}=ℜ{}, i and j are
respectively Dirichlet and Neumann boundary knots. u np is the n-th order of particular solution defined as
∞

∑u

u np =

m
h

.

(16)

m = n +1

The m-order homogeneous solution is approximated by
L

uhm (x ) = ∑ β kmum# (rk ) ,

(17)

k =1

where L is the number of boundary nodes, and u m# is the corresponding m-th order fundamental or
general solutions. Collocating boundary equations (15a,b,c) only on boundary nodes, we have the
boundary discretization equations
L



∑ β u (r ) = R(x ) − u (x ) 
=b
∂u (r )
∂u (x )
(
)
=
−
β
N
x
∑
∂n
∂n 
0
k

k =1

L

0

ik

0
p

i

0

0
p

jk

0
k

i

0

,

(18a)

j

j

k =1


∑ β ℜ {u (r )} = f (x ) − ℜ {u (x )} 
=b
∂ℜ {u (r )} ∂ ( f (x ) − ℜ {u (x )})
β
=
L

1
k

k =1

L

∑

k =1

#
1

n
k

1
p

i

0

jk

i

n
k

1
p

j

,

(18b)

j



∂n

{ }
{ } 
{u (r )} = ∂ (ℜ {f (x )}− ℜ {u (x )}) = b

ℜ n −1 un# (rik ) = ℜ n − 2 { f (xi )} − ℜ n −1 u np (xi )

∂ℜ n −1 n#
β
∑
∂n
k =1
L

0

ik

∂n

k =1

∑β

#
1

1

0

1
k

L

0

n−2

n −1

j

ik

∂n

n
p

n

, n=2,3,….

(18c)

j



In terms of the MRM, the successive process is truncated at some order M, namely, let

{ }

ℜ M −1 u Mp = 0 .

The practical solution procedure is a reversal recursive process:
β kM → β kM −1 → L → β k0 .
It is noted that due to
ℜn −1 {uhn (rk )} = uh0 (rk ) ,
the coefficient matrices of all successive equation are the same, i.e.
Qβ kn = b n , n=M,M-1,…,1,0,
Thus, the LU decomposition algorithm is suitable for this task. Finally, the solution is given by
M

L

u (xi ) = ∑∑ β knun# (rik ) .

(19)
(20)
(21)
(22)
(23)

n = 0 k =1

The BPM can use either singular fundamental solution or nonsingular general solution, respectively
relative to the MFS and BKM. It is noted that the BPM with M=1 degenerates into the BKM or MFS
without using the inner nodes. The only difference between the BKM (MFS) and BPM lies in how to
evaluate the particular solution. The former applies the dual reciprocity principle, while the latter employs
the multiple reciprocity principle. The advantage of the BPM over the BKM is that it dose not require
interior nodes which may be especially attractive in such problems as moving boundary, inverse
problems, and exterior problems. However, the BPM may be more mathematically complicated and
computationally costly due to the iterative use of higher-order fundamental or general solutions. It is
expected that like the MRM[8], the truncated order M in the BPM may not be large (usually two or three

Proceeding of China Congress on Computational Mechanics’ 2001

orders) in a variety of practical uses. The above form of the BPM is unsymmetric. It is rather
straightforward to derive the symmetric BPM scheme[6,7] by replacing Eq. (17) with a representation
similar to Eq. (10) and thus, matrix Q in Eq. (22) will be symmetric if operator ℜ{} is self-adjoint.

3. Numerical experiments
The 2D and 3D irregular geometries tested are illustrated in Figs. 1 and 2, which respectively include a
triangle cut-out and a two-ball (radii=1, center distance= 2 ) cavity. Unless specified Neumann boundary
conditions shown in Fig. 1 and on x=0 surface in Fig. 2, the otherwise boundary are all Dirichlet type.
The BKM employs 9 inner nodes for 2D inhomogeneous case as shown by small crosses in Fig. 1.
Equally spaced knots were applied on the boundary except on two-ball surface where the random knots
were employed. The L2 norms of relative errors are calculated by the numerical solutions at 364 nodes for
2D and 500 nodes for 3D. The absolute error is taken as the relative error if the absolute value of the
solution is less than 0.001. Note that different nodes are used for BKM and BPM coefficients and for L2
norm of relative errors.

Equations of Helmholtz and diffusion-reaction
∇ 2u + γ 2u = f ( x ) ,

and convection-diffusion

∇ 2u − τ 2u = q ( x )

(24a,b)

D∇ 2u ( x ) − v • ∇u ( x) − κu ( x) = g (x ) ,

(25)
are examined, where v denotes a velocity vector, τ is the Thiele parameter, D is the diffusivity coefficient,
κ represents the reaction coefficient. The accurate solutions are
(26)
u = x 2 sin x cos y
for 2D inhomogeneous Helmholtz problem ( γ = 2 ) and
u = sin x cos y cos z
and
u = e −σx + e −σy + e −σz

(27a,b)

for 3D homogeneous Helmholtz ( γ = 3 ) and convection-diffusion (vx=vy=vz=-σ, κ=0, Peclect number is
24 for σ =1 and 480 for σ=20) problems. The corresponding inhomogeneous function f(x), g(x) and
boundary conditions can be derived accordingly. The high-order general solutions of Helmholtz,
diffusion, and convection-diffusion problems are respectively given by
− n 2 +1 + m
− n 2 + 1+ m
um# (r ) = Am (γr )
J n 2 −1+ m (γr ) ,
um# (r ) = Am (τr )
I n 2 −1+ m (τr ) ,
(28)
um# (r ) = Am (µr )

− n 2 +1 + m

v⋅r

e 2 D I n 2 −1+ m (µr ) , n≥2,

(29)

where n is the dimension of the problem; m denotes the order of general solution; J and I respectively
represent the Bessel and modified Bessel function of the first kind; Am=Am-1/(2*m*γ2), A0=1.

[

]

1

µ = ( v 2D) + κ D 2 .
2

(30)

Proceeding of China Congress on Computational Mechanics’ 2001

In this study, the second-order high order general solutions are used as the radial basis function in the
BKM approximation of particular solution. The experimental results are displayed in Tables 1 and 2. It is
found that both methods produce very accurate solutions with a small number of nodes for
inhomogeneous Helmholtz problems. Like the DR-BEM, the BKM is somewhat sensitive to inner node
locations and number. In the 3D case relative to 2D case, more nodes are used due to not only one more
dimension but also larger domain and more complicated boundary shape. It is expected that compared
with other numerical techniques, the BKM and BPM become more efficient for higher-dimensional
complex-shape geometry problems since the general solutions of high-dimensional operators are simpler
and radial basis function is independent of dimensionality and geometry complexity.
Table 1. L2 norm of relative errors for 2D inhomogeneous Helmholtz problems by the BKM and BPM
(numbers inside parentheses indicate boundary plus inner nodes used by the respective method ).
BKM (26+9)
BKM (33+9)
BPM (26)
BPM (33)
1.9e-3
9.3e-5
2.7e-3
6.8e-4
Table 2. L2 norm of relative errors for 3D homogeneous Helmholtz and convection-diffusion problems by
the BKM.
Helmholtz
Convection-diffusion (σ=1)
Convection-diffusion (σ=20)
4.6e-3 (298)
1.7e-4 (466)
9.0e-3 (136)
2.2e-3 (298)
8.8e-15 (136)
6.8e-15 (178)

4. Remarks
The BKM and BPM circumvent the troublesome singular integral inherent in the BEM and are very
easy to learn and program. It is noted that unlike the MR-BEM, the BPM does not need to generate more
than one interpolation matrix, which tremendously reduces the computing effort and storage. In addition,
both schemes are essentially meshfree, spectral convergence and symmetric technique. Similar to the
comparisons between the DR-BEM and MR-BEM, the BKM may be mathematically simpler and more
generally applicable than the BPM, while the latter has advantage not requiring inner nodes for
inhomogeneous problems and is very suitable for problems whose higher-order homogeneous solution
quickly tends to zero. For a complete description and references of the BKM and BPM see refs. [4-7].

References
1. Nardini, D. and Brebbia, C.A., A new approach to free vibration analysis using boundary elements, Applied
Mathematical Modeling, 1983, 7: 157-162.

2. Kansa, E.J. Multiquadrics: A scattered data approximation scheme with applications to computational fluiddynamics. Comput. Math. Appl. 1990, 19: 147-161.

3. Golberg, M.A. and Chen, C.S., The method of fundamental solutions for potential, Helmholtz and diffusion
problems. In Boundary Integral Methods - Numerical and Mathematical Aspects, (Ed. by M.A. Golberg), pp. 103176, Comput. Mech. Publ., 1998.
4. Chen, W. and Tanaka, M. A meshless, exponential convergence, integration-free, and boundary-only RBF
technique. Comput. Math. Appl., (in press), 2001.
5. Chen, W. and Tanaka, M. New Insights into Boundary-only and Domain-type RBF Methods. Int. J. Nonlinear Sci.
& Numer. Simulation, 2000, 1(3): 145-151.
6. Chen, W. Several new domain-type and boundary-type numerical discretization schemes with radial basis function.
ACM Computing Research Repository, http://xxx.lanl.gov/abs/cs.NA/0104018, April, 2001.
7. Chen, W. New RBF collocation schemes and their applications. Int. Workshop Meshfree Methods for Partial
Differential Equations, Bonn, Germany, Sept. 2001.
8. A.J. Nowak and A.C. Neves (ed.). The Multiple Reciprocity Boundary Element Method. Comput. Mech. Publ.,
Southampton, UK, 1994.

