Copyright

c 2005 Tech Science Press CMC, vol.2, no.3, pp.177-188, 2005

The method of fundamental solutions for eigenproblems with Laplace and

biharmonic operators

S.Yu. Reutskiy1

Abstract: In this paper a new meshless method for In particular, the stationary points of the functional
eigenproblems with Laplace and biharmonic operators Z Z
2
in simply and multiply connected domains is presented. R (w) = ∇w dΩ/ w2 dΩ.
The solution of an eigenvalue problem is reduced to a se- Ω

quence of inhomogeneous problems with the differential coincide with eigenfunctions of the Laplace operator.
operator studied. These problems are solved using the See [Courant (1943); Courant and Hilbert (1953); Morse
method of fundamental solutions. The method presented and Feshbach (1953); Strang and Fix (1973)] for more
shows a high precision in simply and multiply connected details and references. Then, using an approximation
domains. The results of the numerical experiments justi- for w with finite number of free parameters, one gets
fying the method are presented. the same problem in a finite-dimensional subspace which
can be solved by a standard procedure of linear alge-
keyword: Method of fundamental solutions, Mem-
bra, e.g., see [Golub and Loan (1996); Strang (1976)].
branes and Plates, Free vibration problem
However, a standard finite differences method can pro-
duce good results when dealing with a particular type
1 Introduction of shapes defined on rectangular grids, while for other
The goal of this paper is to present a new numerical tech- type of shapes the finite element method or the bound-
nique for solution of the following eigenvalue problems: ary element method are more appropriated. The method
of fundamental solutions (MFS) [Fairweather and Kara-
∇2 w + k2 w = 0, x ∈ Ω ⊂ R 2 , (1) georghis (1998); Golberg and Chen (1998, 1997)] is the
convenient tool in this field. The similar technique is
B1 [w] = 0, x ∈ ∂Ω
used in the boundary knot method (BKM) [Chen (2005);
and Chen and Tanaka (2002)]. Unlike the MFS, it employs
nonsingular general solutions as the basis functions to
∇4 w − k4 w = 0, x ∈ Ω ⊂ R 2 , (2) avoid the fictitious boundary outside the physical do-
main.
w = 0, B2 [w] = 0, x ∈ ∂Ω.
In the framework of the boundary methods a general ap-
Here Ω is a simply or multiply connected domain of in- proach to solving these problems is as follows. First,
terest with boundary ∂Ω. The boundary operator in (1) using an integral representation of w in the BEM, or an
B1 [...] will be considered of the two types: the Dirichlet approximation over fundamental solutions in MFS, one
B1 [w] = w and of the Neumann type B1 [w] = ∂w/∂n; for gets a homogeneous linear system A (k) q = 0 with ma-
biharmonic operator in (2), B2 [w] = ∂w/∂n or B2 [w] = trix elements depending on the wave number k. The de-
∂2 w/∂n2 . As a mechanical application, this corresponds terminant of this matrix must be zero to obtain the non-
to recovering the free vibration frequencies of mem- trivial solution:
branes and plates. Such problems often arise in engi-
neering applications. det[A (k)] = 0 (3)
The usual approach for eigenvalue problems with a self- This equation must be investigated analytically or numer-
adjoint operator is to use the Rayleigh minimal principle.
ically to get the eigenvalues. This technique is described
1 Laboratory of Magnetohydrodynamics, Timurovtzev, 29 D, ap. in [Karageorghis (2001); Chen, Lin, Kuo, and Chyuan

51, 61142, Kharkov, Ukraine (2001); Chen, Liu, and Hong (2003); Chen, Chen, Chen,
178 Copyright

c 2005 Tech Science Press CMC, vol.2, no.3, pp.177-188, 2005

Lee, and Yeh (2004); Chen, Chen, and Lee (2005)] with following 1D Sturm − Liouville problem on the interval
more details. In the two latest papers there is a complete [0, 1]:
bibliography on the subject considered. In [Alves and
2
Antunes (2005)] the MFS based procedure is applied to d w + k2 w = 0, w (0) = w (1) = 0. (6)
eigenproblems in general simply connected shapes. dx2
The method presented in this article uses the same MFS The well known solution is: kn = nπ, wn = sin(nπx),
boundary technique. This is a mathematical model of n = 1, 2, ..., ∞.
physical measurements when the resonance frequencies Following the boundary approach, let us consider the
of a system are determined by the amplitude of response fundamental solution
to some external excitation. As a result, e.g., instead
of (1) we solve a sequence of inhomogeneous boundary 1
Ψ (x, ξ, k) = exp (ik |x − ξ|), (7)
value problems (BVP): 2k
which satisfies the homogeneous equation everywhere
∇2 w + k2 w = f (x), x ∈ Ω ⊂ R 2 , (4)
except the singular point x = ξ. A general solution of
B1 [w] = 0, x ∈ ∂Ω, the homogeneous equation in the interval [0, 1] can be
written in the form:
where f describes some source placed outside the solu-
tion domain. Let F (k) be some norm of the solution w. w = q1 Ψ (x, ξ1 , k) + q2 Ψ (x, ξ2 , k) .
As it will be shown below, this function of k has sharp
maximums at the eigenvalues and, under some condi- Here ξ1 , ξ2 are two source points placed outside the solu-
tions can be used for their determining. Certainly such tion domain [0, 1], e.g., ξ1 < 0, ξ2 > 1; q1 , q2 are free pa-
behaviour of F (k) near the eigenvalues is a consequence rameters. Using the boundary conditions w (0) = w (1) =
of (3). Techniques of numerical solution of linear BVPs 0, one gets the linear system:
like (4) are well developed. It should be emphasized that

any Helmholtz (or biharmonic) equation solver can be q1 e(−ikξ1 ) + q2 e(ikξ2 ) = 0
A (k) q =
used in the framework of the method presented. How- q1 e(−ik(1−ξ1 )) + q2 e(ik(ξ2 −1)) = 0
ever, the MFS technique seems to be a more suitable one
for this goal in the case of an arbitrary domain. The wave numbers kn can be determined from the con-
The outline of this paper is as follows: for the sake of dition: det[A (k)] = 0. After simple transforms we get:
simplicity we begin by describing the 1D case in Section exp (2ik) = 1, or k = nπ. Thus, MFS gets the exact solu-
2. In Section 3, we present the algorithm of MFS in ap- tion. Note that in multidimensional cases such computa-
plication to problem (1). Here we present numerical ex- tions are time consuming and not so simple.
amples to illustrate the method presented for simply and As it is mentioned above, the method suggested is a
multiple connected domains. In Section 4, the same tech- mathematical model of physical measurements, when a
nique is described in application to problem (2). Some mechanical or acoustic system is excited by an external
generalization of the technique and the fields of its de- source and resonance frequencies can be determined us-
velopment are discussed in Section 5. ing an increase of amplitude of oscillations near these
frequencies. So, instead of (6) we solve the inhomoge-
2 One-dimensional eigenproblem neous problem:

To illustrate the method presented let us consider the d2w

+ k2 w = f (x), w (0) = w (1) = 0. (8)
wave equation [Morse and Feshbach (1953)] dx2
∂2 u ∂2 u The general solution can be written in the form:
= 2 (5)
∂t 2 ∂x
w = q1 Ψ (x, ξ1 , k) + q2 Ψ (x, ξ2 , k) + w p . (9)
with the Dirichlet conditions at the endpoints of the in-
terval [0, 1], i.e., u (0,t) = u (1,t) = 0. Considering the When the excitation is performed by the point source
free harmonic vibrations u (x,t) = e−ikt w (x), we get the with the same wave number k which is placed at the point
The method of fundamental solutions for eigenproblems with Laplace and biharmonic operators 179

ξ0 outside the solution domain, then f (x) = iδ (x − ξ0 )

lnFd
and the particular solution is:
1 15
w p = Ψ (x, ξ0 , k) = exp (ik |x − ξ0 |) . (10)
2k
10
Using again the same homogeneous boundary conditions
w (0) = w (1) = 0, now we get an inhomogeneous linear
5
system for each k. Let us introduce the norm of the solu-
tion as
 k
1 N 5 10 15 20 25 30
F (k) = ∑ |w (xn)|2, Fd (k) = F(k)/F(k0), (11)
N n=1
Fd
where Fd (k) is the dimensionless value, k0 is a reference 10
wave number and the points xn are randomly distributed 8
in [0, 1]. In all the calculations presented in this section
6
we use N = 5. This function characterizes the value of
the response of the system to the outer excitation. Note 4
that the right hand side f corresponding to (10) equals
2
to zero identically inside [0, 1] and BVP (8) has a unique
solution w = 0 for all k except k = kn - eigenvalues when k
3 4 5 6 7
the solution is not unique.
In Fig. 1 the value Fd as a function of the wave number k Figure 1 : Resonance curve in 1D eigenproblem.
is shown.
The graph contains large sharp peaks at the positions of
eigenvalues. Generally speaking, this resonance curve lnFd
can be used to determine the eigenvalues in the same way
as det [A (k)] in the technique described above. However,
the graph Fd (k) is a non smooth one, as it is shown in 15
the lower part of Fig. 1 with more details. This can be
10
explained by the following reasons. Problem (8), (9)
with w p given in (10) has the exact solution q1 = 0, 5
q2 = −eik(ξ0 −ξ2 ) and so the total solution w(x) = 0, for x ∈
[0, 1]. So, here we have F (k) which is equal to zero with k
machine precision accuracy when k is far from eigenval- 5 10 15 20 25 30
ues; F (k) grows considerably in a neighbourhood of the
eigenvalues when the linear system becomes almost de- Fd
generated. And a smoothing procedure is needed to get
an appropriate curve which is convenient for applying an 15
optimization procedure. The following two smoothing
procedures are used in the paper. 10

2.1 smoothing by a dissipative term 5

The first procedure consists of introducing an additional
k
dissipative term in the governing equation. And instead 3 4 5 6 7
of (8) we consider the problem:
Figure 2 : Resonance curve in 1D eigenproblem.
d2w  2 
2
+ k + iεk w = f , w (0) = w (1) = 0. (12) Smoothing by a friction term.
dx
180 Copyright

c 2005 Tech Science Press CMC, vol.2, no.3, pp.177-188, 2005

Here ε is a small parameter. This means that the ini- Table 1 : One dimensional eigenproblem. The relative
tial wave equation (5) is changed by the equation ∂tt2 u = errors in calculations of the eigenvalues. Smoothing by
∂2xx u − ε∂t u which describes vibration of a homogeneous the friction
(ex)
term.
string with a friction [Morse and Feshbach (1953)]. The ki ε = 0.1 ε = 10−3 ε = 10−6
fundamental solution is: π 1.3 · 10−4 1.3 · 10−8 1.7 · 10−12
2π 3.2 · 10−5 3.1 · 10−9 1.6 · 10−12
1 3π 1.4 · 10−5 1.4 · 10−9 1.5 · 10−12
Ψ (x, ξ, k, ε) = exp (iχ |x − ξ|) ,
2χ 4π 7.9 · 10−6 7.9 · 10−10 9.7 · 10−13

χ = k2 + iεk. (13) 5π 5.1 · 10−6 5.0 · 10−10 9.0 · 10−13

Now the system w (0) = 0, w (1) = 0 with w p given in

(10) has a unique non zero solution for all real k. The
resonance curve corresponding to ε = 10−6 is shown in 2.2 smoothing by a shift of wave numbers
Fig. 2.
The second smoothing technique is as following. Let
Now this is a smooth curve with separated maximums
us introduce the constant shift ∆k between the exciting
at the positions of eigenvalues. To find the eigenvalues
source and the studied mode, i.e., instead of (10), we take
we use the following algorithm through the paper. Let us
the particular solution in the form:
look for the eigenvalues on the interval [a, b] Then:
(A) w p = Ψ (x, ξ0 , k + ∆k)
step 0: Choose h > 0; 1
if F (a) > F (a + h) goto step 5; = exp (i (k + ∆k) |x − ξ0 |) . (15)
2 (k + ∆k)
step 1: x1 = a; F1 = F (x1 ) ;
step 2: x2 = x1 + h; F2 = F (x2 ) ; Now the linear system w (0) = w (1) = 0 has non zero so-
if x2 > b stop; lutions for all k except the eigenvalues kn when the sys-
step 3: if F2 > F1 then [F1 = F2; x1 = x2 ]; tem becomes degenerate. However, due to iterative pro-
goto step 2; cedure of solution and rounding errors we never solve the
step 4: find the maximum point xm of F (x) system with the exact kn . And we observe degeneration
on [x2 − 2h, x2 ] ; of the system as a considerable growth of the solution in a
step 5: x1 = a; F1 = F (x1 ) ; neighbourhood of the eigenvalues. The resonance curve
step 6: x2 = x1 + h; F2 = F (x2 ) ; corresponding to ∆k = 1 is shown in Fig. 3.
if x2 > b stop; Some results of the calculations we got using the second
step 7: if F2 < F1 then [F1 = F2; x1 = x2 ; smoothing technique are presented in Tab. 2. The values
goto step 6]; ξ1 , ξ2 , ξ0 are the same as above.
else goto step 2.
Below we will name these procedures as ε−procedure
Note that any univariate optimization procedure can be and k−procedure.
used at step 4. In particular, we applied Brent’s method
based on a combination of parabolic interpolation and bi-
section of the function near to the extremum(see [Press, Table 2 : One dimensional eigenproblem. The relative
Teukolsky, Vetterling, and Flannery (2002)], Ch. 10 and errors in calculation of the eigenvalues. Smoothing by
[Brent (1973)], Ch. 5 ). The step is taken h = 0.01 shift of the wave numbers.
through the paper. The data placed in Tab. 1 are ob- (ex)
ki △k = 0.1 △k = 1 △k = 10
tained by applying this technique with ε = 0.1, 10−3 , π 1.4 · 10−11 9.1 · 10−12 7.8 · 10−12
10−6 . Other parameters are: ξ1 = −0.5, ξ2 = 1.5, ξ0 = 5. 2π 5.8 · 10−13 3.5 · 10−12 5.5 · 10−12
Here we place the relative errors 3π 6.4 · 10−12 1.3 · 10−12 3.5 · 10−12
(ex) (ex) 4π 3.3 · 10−13 2.8 · 10−12 2.3 · 10−12
er = |ki − ki |/ki (14)
5π 5.3 · 10−12 3.5 · 10−12 5.9 · 10−13
in the calculation of the first five eigenvalues.
The method of fundamental solutions for eigenproblems with Laplace and biharmonic operators 181

Here the points xi , i = 1, ..., Nc are distributed uniformly

lnFd
on the boundary. We take Nc approximately twice as
large as the number of free parameters N. The problem
30 is solved by the standard least squares procedure. Note
that we get (18) as a result of discretization of the integral
20 condition:
Z
min {B1 [w (x|q)]}2 ds
10 w ∂Ω
More details of this technique can be found, e.g., in [Fair-
k weather and Karageorghis (1998); Golberg and Chen
5 10 15 20 25 30
(1998)].
Fd As a particular solution corresponding to the exciting
source we take the same fundamental solution
(1)
15 w p (x) = Φex(x, ζex , k) ≡ H0 (k |x − ζex |) (19)
with ζex placed outside the solution domain.
10
When dealing with problems in multiply connected do-
mains, the same trial functions can be used. And the
5
source points should be placed also inside each hole. As
k
an alternative approach one can use the special trial func-
3 4 5 6 7 tions associated with each hole:
(1)
Φs,1 (x) = H0 (krs),
Figure 3 : Resonance curve in 1D eigenproblem. (1)
Φs,2n+1 (x) = Hn (krs) cosnθs ,
Smoothing by a shift between the wave numbers. (1)
Φs,2n (x) = Hn (krs ) sinnθs . (20)

Here rs = |x − xs |, θs is the local polar coordinate sys-

3 Helmholtz eigenproblem tem with the origin at xs . This is so-called Vekua ba-
Applying the MFS to problem (4) we look for an approx- sis [Vekua (1957); Hafner (1990)], or multipole expan-
imation solution in the form of a linear combination: sion. It is proven that every regular solution of the 2D
N Helmholtz equation in a domain with holes can be ap-
w (x|q) = w p (x) + ∑ qn Φn (x), (16) proximated with any desired accuracy by linear combi-
n=1 nations of such functions if the origin xs of a multipole is
where w p is the particular solution corresponding to f , inside every hole. In this case instead of (16) we use:
and the trial functions N
(1) w (x|q, ps ) = w p (x) + ∑ qn Φn (x)
Φn (x) = H0 (k |x − ζn |) (17)
n=1
satisfy the homogeneous PDE. This is the so-called S M

Kupradze basis [Kupradze (1967)]. The singular points +∑ ∑ ps,m Ψs,m (x), (21)
s=1 m=1
ζn are located outside the solution domain. The free
parameters qn should be chosen to satisfy the boundary where S is the number of holes and M is the number of
condition B1 [w (x|q)] = 0, x ∈ ∂Ω. In particular the un- terms in each multipole expansion.
knowns qn are taken as a solution of the minimization When the ε−smoothing procedure is applied, then in-
problem: stead of (4) we consider the problem:
 2  
Nc N ∇2 wh + k2 + iεk wh = 0, x ∈ Ω,
min ∑ B1 [w p (xi )] + ∑ qn B1 [Φn (xi )] (18)
q
i=1 n=1 B1 [wh (x)] = −B1 [w p (x)], x ∈ ∂Ω. (22)
182 Copyright

c 2005 Tech Science Press CMC, vol.2, no.3, pp.177-188, 2005

with some small ε > 0. Note that this problem has a Table 3 : Circular domain with the radius r = 1.
unique nonzero solution for all real k. Then the trial func- The relative errors in calculations of the eigenvalues.
tions (17) should be also modified: ε−procedure; ε = 10−6 .
Dirichlet condition
(1)
Φn (x) = H0 (χ |x − ζn |), i N = 15 N = 20 N = 25
 1 8 · 10 −11
8 · 10 −12
7 · 10−12
χ (k, ε) = k2 + iεk. (23) −3 −11
2 2 · 10 5 · 10 2 · 10−11
3 3 · 10−9 1 · 10−9 1 · 10−9
Applying the k−procedure we modify the particular so- −3
4 2 · 10 4 · 10−11 1 · 10−11
lution which should be taken in the form:
5 6 · 10−7 2 · 10−3 1 · 10−9
(1)
w p (x) = Φex (x, k) ≡ H0 (k |x − ζex |), 
k = k + ∆k. (24) Neumann condition
1 2 · 10−9 2 · 10−9 2 · 10−9
−9
3.1 numerical examples 2 4 · 10 2 · 10−9 2 · 10−9
3 9 · 10−12 1 · 10−11 6 · 10−12
Here the results of the numerical experiments are given
4 7 · 10−8 9 · 10−10 8 · 10−10
to illustrate the method presented. In all the cases consid-
5 2 · 10−6 6 · 10−10 3 · 10−10
ered below the resonance curve F (k) is computed using
Nt testing points xt,l ∈ Ω: F (k) = 1/Nt ∑Nl=1 t
|w (xt,l )|2 . Table 4 : Circular domain with the radius r = 1. Dirich-
In all the calculations we use 15 testing points distributed let condition. The relative errors in calculations of the
inside Ω with the help of RNUF generator of pseudoran- eigenvalues. ε−procedure with varying ε, N = 25.
dom numbers from the Microsoft IMSL Library. To get
the eigenvalues we look for the maxima of F (k) using
the Brent’s procedure mentioned. i ε = 10−2 ε = 10−4 ε = 10−6 ε = 10−8
Example 1) A circular domain with the radius r = 1 sub- 1 6.4 · 10−6 6.0 · 10−10 7.3 · 10−12 4.9 · 10−11
jected to Dirichlet or Neumann boundary condition is 2 2.4 · 10−6 1.9 · 10−10 2.0 · 10−11 4.3 · 10−11
considered. The exciting source is placed at the posi- 3 3.2 · 10−6 1.4 · 10−9 1.0 · 10−9 −
4 9.0 · 10 −7 1.6 · 10 −10 1.3 · 10 −11 −
tion ζex = (5, 5); the singular points ζn of the fundamen-
−6 −9 −9
tal solutions (17) are located on the circle with the ra- 5 1.1 · 10 1.6 · 10 1.4 · 10 −
dius R = 2. The results shown in Tab. 3 correspond to 6 6.5 · 10−7 1.5 · 10−10 − −
−6 −7 −10
ε = 10 . Here we place the relative errors (14) in the 7 4.9 · 10 4.8 · 10 − −
calculation of the first 5 eigenvalues. The line − in a cell 8 2.7 · 10−6 1.1 · 10−9 − −
indicates that the solution process failed with these pa- 9 4.9 · 10 −7 5.9 · 10 −9 − −
(ex)
rameters. The exact eigenvalues ki are the roots of the 10 5.2 · 10−6 − − −
′
equation Jn (k) = 0 (Dirichlet) or Jn (k) = 0 (Neumann).
Example 2) The role of the parameter ε is shown in
Tab. 4. We solve the same problem as above with Dirich-
let condition. Here we fix the number of free parame- sharply part of the resonance curve. As a result the al-
ters N = 25 and vary the parameter ε. The parameter ε gorithm (A) ’jumps over’ the eigenvalues and one should
coarsens the system. For a large ε we can calculate all thedecrease the step parameter h to capture the maxima. As
it is shown in Tab. 4, for ε = 10−8 the algorithm finds k1
eigenvalues ki , i = 1, ..., 10 but the precision is not very
high. When ε decreases, the precision in determining of and k2 with h = 0.01. When h is reduced to 0.001 then
ki increases but it fails for large i. the algorithm also gives the eigenvalues k3 and k4 . To get
The figures Fig. 4, Fig. 5, Fig. 6 correspond to the data ki , i = 1, ..., 10 one should take h = 0.0001. However, the
placed in Tab. 4. For ε = 10−2 the resonance peaks are algorithm becomes highly expansive in the CPU time.
spread because the friction. When ε decreases the peaks Example 3) Next, we consider the case when Ω is the
become more sharp and narrow. Besides for ε = 10−8 the unit square with the same Dirichlet or Neumann bound-
peaks corresponding to ki , i > 2 are placed on the rising ary condition.This problem has an analytical solution:
The method of fundamental solutions for eigenproblems with Laplace and biharmonic operators 183

Fd the results of calculation of the first 5 eigenvalues with

500 ε = 10−6 . The placement of the singular points ζn and
the exciting source are the same as above.
400
Example 4) For the next example, we consider an an-
300 nular case of the doubly connected domain between the
two circles: Ω = (x1 , x2 ) | r12 ≤ x21 + x22 ≤ r22 The in-
200 ner and outer radii of an annular domain are r1 = 0.5
100 and r2 = 2, respectively. We take Dirichlet condition
on the outer boundary and Neumann on the inner one.
k The singular points are distributed at the circles with
3 4 5 6
the radii a = 5(outside the domain) and b = 0.3 (inside
the hole). The number of the singular points on each
Figure 4 : Circular membrane with the radius 1. Dirich- auxiliary contour is equal to N. The exciting source is
let conditions. ε - procedure with ε = 10−2 . placed at ζex = (10, 10). In Tab. 6 we present the rela-
tive errors (14) in calculation of the first 5 eigenvalues of
Fd the problem described with ε = 10−5 . The values ki
(ex)
500
are obtained numerically as the roots of the equation:
400 Jn′ (r1 k)Yn (r2 k) − Jn (r2 k)Yn′ (r1 k) = 0.

300
Table 5 : Square with the side a = 1. The relative er-
200 rors in calculations of the eigenvalues. ε−procedure;
100 ε = 10−6 .
Dirichlet condition
k i N = 15 N = 20 N = 25
3 4 5 6
1 1 · 10−6 3 · 10−8 1 · 10−9
2 1 · 10−5 9 · 10−8 1 · 10−8
Figure 5 : Circular membrane with the radius 1. Dirich- 3 8 · 10−5 3 · 10−8 8 · 10−9
let conditions. ε - procedure with ε = 10−4 . 4 3 · 10−4 1 · 10−6 3 · 10−9
5 3 · 10−3 4 · 10−5 6 · 10−7
Fd
500 Neumann condition
1 4 · 10−7 5 · 10−8 8 · 10−12
400 2 1 · 10−6 3 · 10−8 3 · 10−9
300 3 4 · 10−5 1 · 10−7 3 · 10−10
4 1 · 10−4 6 · 10−6 5 · 10−9
200 5 5 · 10−4 2 · 10−5 6 · 10−7
100 Table 6 : Annular domain. The relative errors in calcu-
k lations of the eigenvalues. ε−procedure; ε = 10−5 .
3 4 5 6
(ex)
i ki N = 15 N = 20 N = 25
Figure 6 : Circular membrane with the radius 1. Dirich- 1 1.3339427880 5 · 10−11 2 · 10−11 2 · 10−11
let conditions. ε - procedure with ε = 10−8 . 2 1.7388632616 6 · 10−8 7 · 10−12 5 · 10−12
3 2.4753931967 − 7 · 10−11 8 · 10−12
4 3.1645013237 − 7 · 10−8 5 · 10−11

k(ex) = π i2 + j 2 , i, j = 1, 2, ..(Dirichlet condition) and 5 3.2899912986 − − 7 · 10−11
i, j = 0, 1, 2, ..(Neumann condition). In Tab. 5, we show
184 Copyright

c 2005 Tech Science Press CMC, vol.2, no.3, pp.177-188, 2005

Example 5) In this example, doubly connected region Table 7 : Circle with a small hole. Dirichlet boundary
with the inner region of vanishing maximal dimension is condition. The outer radius: r2 = 2. The relative errors
considered. The geometry of the problem is the same in calculation of the first ten eigenvalues. k−procedure
as in Example 3. However, here we consider the case with ∆k = 1.
r1 = 0.1, N = 50, M = 11
of very small inner holes. In particular, we take r1 = (ex)
10−1 , 10−2 , 10−3 with the same fixed r2 = 2. Now, the i ki er
Kupradze type basis functions (17) are unfit to approxi- 1 1.5322036536 1.9 · 10−8
mate the solution in a neighbourhood of the hole. Here 2 1.9301625755 5.8 · 10−9
we use a combined basis which includes the trial func- 3 2.5680354360 1.6 · 10−9
tions (17) with the singular points placed on an auxiliary 4 3.1900833197 1.3 · 10−11
circular contour outside the solution domain and a mul- 5 3.2126996563 7.4 · 10−9
tipole expansion with the origin at the center of the hole. 6 3.5522743165 3.7 · 10−10
Thus, we look for an approximate solution in the form: 7 3.7941712382 1.2 · 10−11
N M 8 4.2101115868 9.0 · 10−12
w (x|q, p) = w p (x) + ∑ qn Φn (x) + ∑ pm Ψm (x). 9 4.3857419081 4.4 · 10−12
n=1 m=1 10 4.8805392651 1.0 · 10−11
The data presented in Tab. 7, Tab. 8, Tab. 9 correspond
to the number of sources on the outer auxiliary circu- Table 8 : Circle with a small hole. Dirichlet boundary
lar contour N = 50. The number of terms in multi- condition. The outer radius: r2 = 2. The relative errors
pole expansion M varies from M = 11(r1 = 10−1 ) to in calculation of the first ten eigenvalues. k−procedure
M = 5(r1 = 10−3 ). The exciting source is placed at with ∆k = 1.
the position ζex = (10, 10). We use the k−procedure
r1 = 0.01, N = 50, M = 7
with the shift ∆k = 1. We would like to draw the read- (ex)
ers’ attention to the fact that the method presented can i ki er
(ex)
separate very close eigenvalues: k4 = 3.1900833197 1 1.3709447159 2.5 · 10−8
(ex) 2 1.9160005377 5.4 · 10−9
and k5 = 3.2126996563(see data corresponding to r1 =
3 2.5678112121 1.6 · 10−9
10−1 ). Here the step in the algorithm (A) is taken
h = 0.001. The detailed discussion of Vecua basis for 4 2.9632630840 5.3 · 10−9
Helmholtz equation can be found in [Hafner (1990)]. 5 3.1900809955 2.9 · 10−12
6 3.5082790790 2.3 · 10−12
4 Eigenproblems with biharmonic operator 7 3.7941712738 1.0 · 10−9
8 4.2086222910 7.6 · 10−12
According to the technique proposed, instead of (2) let’s 9 4.3857419733 1.1 · 10−11
consider BVP 10 4.5543927267 1.3 · 10−9
∇4 w − k4 w = f , x ∈ Ω ⊂ R 2 , (25)
w = 0, B2 [w] = 0, x ∈ ∂Ω.
(1)
In application to this problem, the MFS technique is sim- where H0 is the Hankel function and K0 is the modi-
ilar to the one considered in the previous section. The fied Bessel function of the second kind and of order zero.
trial functions now are of the two types: the fundamental So, an approximate solution is sought in the form of the
solutions of the Helmholtz operator ∇2 + k2 : linear combination:
(1) (1) N N
Φn (x) = H0 (k |x − ζn |) (26) w (x|q , q ) = w (x) + (1) (2)
1 2 p ∑ q1,nΦn (x) + ∑ q2,nΦn (x).
n=1 n=1
considered above and the fundamental solutions of the
(28)
modified Helmholtz operator ∇2 − k2 :
2 where w p (x) is a particular solution corresponding to the
(2) (1)
Φn (x) = H0 (ik |x − ζn |) = −i K0 (k |x − ζn |) , (27) external source. We take it in the same form as in the
π
The method of fundamental solutions for eigenproblems with Laplace and biharmonic operators 185

Table 9 : Circle with a small ole. Dirichlet boundary should be replaced by the following one:
condition. The outer radius: r2 = 2. The relative errors
 
in calculation of the first ten eigenvalues. k−procedure ∇4 w (x) − k4 + iεk2 w (x) = f (29)
with ∆k = 1.
r1 = 0.001, N = 50, M = 5 (1)
(ex)
and so the arguments of the trial functions Φn (x),
i ki er (2)
Φn (x) should be modified. Applying the k−procedure
1 1.3148533741 2.0 · 10−8 we modify the external source and take it in the form
2 1.9158544900 5.4 · 10−9 (24).
3 2.5678111892 1.5 · 10−9
4 2.8883437835 2.8 · 10−9 4.1 numerical examples
5 3.1900809955 1.1 · 10−10
Example 6) A circular plate with the radius r = 1 sub-
6 3.5077982552 3.0 · 10−11
jected to the boundary conditions: a) w = ∂w/∂n = 0
7 3.7941712738 1.2 · 10−11
(clamped boundary) and b) w = ∂2 w/∂n2 = 0 is con-
8 4.2086221329 5.9 · 10−12
sidered. The exciting source is placed at the position
9 4.3857419733 1.2 · 10−12
ζex = (5, 5); the singular points ζn of the fundamental so-
10 4.4650868082 3.6 · 10−10
lutions (26), (27) are located on the circle with the radius
Table 10 : A circular plate with the radius: r = 1;the R = 2. Remark that now the number of free parameters
relative errors in calculation of the first eight eigenvalues. is 2N. The data presented in Tab. 10 are obtained using
k−procedure, ∆k = 0.1. k−procedure with ∆k = 0.1. Here we place the relative
(ex)
errors (14). The exact eigenvalues ki are the roots of
′ ′
w = ∂w/∂n = 0 the equation Jn (k)In (k) − Jn (k) In (k) = 0 (conditions a))
i N = 20 N = 25 N = 30 or Jn′′ (k) In (k) − Jn (k) In′′ (k) = 0 (conditions b)).
1 3 · 10−9 3 · 10−9 3 · 10−9 Example 7) Next, we consider a square plate with the
−7 −9
2 7 · 10 3 · 10 2 · 10−9 side a = 1 subjected to the boundary conditions w =
−5 −8 −10 ∂2 w/∂n2 =
3 1 · 10 8 · 10 6 · 10  0. This problem has an analytical solution:
4 6 · 10 −10
6 · 10 −10
6 · 10 −10 k(ex) = π i2 + j2 , i, j = 1, 2, .. The results placed in
5 7 · 10−5 9 · 10−7 1 · 10−8 Tab. 11 are obtained using k−procedure with ∆k = 0.1.
6 1.2 · 10 −5 1 · 10 −7 2 · 10 −9
Example 8) A rectangular 1.2 × 0.9 plate subjected to the
7 9.4 · 10 −4 8 · 10 −6 8 · 10−8 boundary conditions w = ∂w/∂n = 0 (clamped bound-
8 1.3 · 10−4 2 · 10−6 2 · 10−8 ary) is considered. The results placed in Tab. 12 are ob-
2 2
w = ∂ w/∂n = 0 tained using k−procedure with ∆k = 0.1. In this case, the
1 2 · 10−9 3 · 10−9 3 · 10−9 analytic solution is not available. The results obtained
2 5 · 10−5 3 · 10−6 1 · 10−7 in [Chen, Chen, Chen, Lee, and Yeh (2004)] and [Kang
−4 −5
3 3 · 10 2 · 10 8 · 10−7 and Lee (2001)] are used for comparison. These data are
−9 −9 −9 placed in the last two columns of the table. Note that us-
4 4 · 10 5 · 10 4 · 10
5 2 · 10 −3
9 · 10 −5
4 · 10 −6 ing ε−procedure with ε = 0.01 and N = 56, we get the
6 3 · 10 −4 1 · 10 −5 6 · 10 −7 following eigenvalues: k1 = 5.95263, k2 = 7.70983, k3 =
7 9 · 10 −3 3 · 10 −4 2 · 10 −5 9.12854, k4 = 10.27133, k5 = 11.96763, k6 = 12.49617
8 9 · 10 −2 6 · 10 −5 3 · 10 −6

5 Concluding remarks

In this paper, a new meshfree method for eigenproblems

with Laplace and biharmonic operators is proposed. This
previous section, i.e. (19). The free parameters are deter- is a mathematical model of physical measurements, when
mined from the boundary conditions. a mechanical or acoustic system is excited by an exter-
We apply the same ε and k smoothing procedure. nal source and resonance frequencies can be determined
When the ε−procedure is applied the governing equation using the growth of amplitude of oscillations near these
186 Copyright

c 2005 Tech Science Press CMC, vol.2, no.3, pp.177-188, 2005

Table 11 : A square plate. The relative errors in calcula- the solution domain.
tion of the first six eigenvalues. k−procedure, ∆k = 0.1.
To test BKM in the framework of the method presented
i N = 20 N = 25 N = 30 N = 40
we solve the same problem as the one described in Ex-
1 3.3 · 10−6 2.0 · 10−7 1.8 · 10−8 1.7 · 10−8
ample 1 with Dirichlet condition. The half of the source
2 6.9 · 10−4 1.5 · 10−5 3.2 · 10−8 1.7 · 10−8
points ζn , n = 1, ..., 1/2N are placed uniformly on the
3 − 7.9 · 10−5 3.7 · 10−6 1.7 · 10−8
boundary ∂Ω. The rest source points ζn , n = 1/2N +
4 − 9.2 · 10−5 5.5 · 10−7 1.7 · 10−81, ..., N are distributed inside Ω with the help of the gen-
5 − 4.5 · 10−2 1.5 · 10−5 2.2 · 10−8erator of pseudorandom numbers. The data presented in
6 − − 1.4 · 10−3 1.2 · 10−7
Tab. 13 are obtained using k−procedure with ∆k = 0.1.
The parameters of the exciting source are the same as
Table 12 : A rectangular plate 1.2 × 0.9 with clamped
above in Example 1.
boundary.
It should be noted that the BKM and the MFS, as well as
i N = 35 N = 42 N = 49 I II the all methods of the Trefftz type in general, have a nar-
1 5.9515 5.9529 5.9527 5.952 5.952 row field of application. It is restricted by the cases when
2 7.7125 7.7116 7.7104 7.703 7.703 there exists a representative set of known exact solutions
3 9.1333 9.1319 9.1316 9.129 9.131 of PDEs under consideration, i.e. by the problems posed
4 9.9466 9.9510 9.9493 9.947 9.955 by linear PDEs with constant coefficients. See, however,
5 10.2692 10.2717 10.2742 10.266 10.27 [Reutskiy (2002)], where a Trefftz type technique is de-
6 11.9501 11.9552 11.9565 11.95 11.95 veloped for PDEs with varying coefficients.
7 12.3849 12.3719 12.3710 − − Besides the Trefftz type techniques produce the systems
of equations with unsymmetric fully populated matri-
Table 13 : The BKM solution. Circular domain with ces. As a result, the MFS is highly ill conditioned. In
Dirichlet conditions. The relative errors in calculations some cases one can overcomes this drawback by the use
of the eigenvalues. k−procedure; ∆k = 0.1. of matrices of the special block circulant structure and
an efficient matrix decomposition technique [Tsangaris,
i N = 10 N = 14 N = 20 N = 30 Smyrlis, and Karageorghis (2004)].
1 2 · 10−4 2 · 10−6 4 · 10−9 7 · 10−9 However, taking in mind further applications of the
2 3 · 10−4 4 · 10−7 1 · 10−10 1 · 10−8 method presented in the paper to eigenproblems with
3 − 9 · 10−5 2 · 10−8 1 · 10−8 PDEs of general type in irregular domains, one should
4 − − 4 · 10−7 4 · 10−9 combine it with meshless methods based on the local
5 − − 1 · 10−6 8 · 10−9 approximation of the solution like the Meshless Local
Petrov-Galerkin Method [Atluri (2004), Han and Atluri
(2003), Han and Atluri (2004)]. The comparison be-
tween global and local approximation, e.g. BEM and
frequencies. The method shows a high precision in sim- FEM, and they combination see in [Grannell and Atluri
ply and multiply connected domains. The idea can be (1978)].
extended quite simply to the 3D case.
Comparing the method with the technique based on com-
The method presented is based on the MFS solution of putations of the determinant of the system, the following
the problem. However, it can be combined with other circumstances should be taken into account. Since the
boundary techniques. The BKM mentioned in Section 1 MFS is highly ill conditioned, the determinant is very
seems to be perspective in this connection. For example, small. Indeed, let us consider again the same eigenvalue
if the BKM is applied to Helmholtz equation, the approx- problem which is described in Example 1 , i.e. Helmholtz
imation solution is looked for in the form: equation in the circle with the radius 1 and Dirichlet
N boundary condition. We take the number of the sources
w (x|q) = w p (x) + ∑ qn J0 (k |x − ζn |) N equal to the number of the collocation points on the
n=1
boundary. Thus, we get a square matrix of the problem
cf. (16). Here the source points ζn can be placed inside
The method of fundamental solutions for eigenproblems with Laplace and biharmonic operators 187

Table 14 : Circular domain with Dirichlet conditions. The number of the source points N = 30; ε−procedure.

ε = 10−1 ε = 10−4 ε = 10−6

i er F(ki ) er F(ki ) er F(ki )
1 4×10−4 0.701 4×10−10 0.701 5×10−12 0.701
2 2×10−4 0.652 1×10−10 0.654 6×10−11 0.654
3 9×10−5 0.509 9×10−10 0.516 1×10−9 0.516

A(k, N) and can calculate the determinant |detA (k, N)| . Brent, R. P. (1973): Algorithms for minimization with-
Placing the sources on the circle with radius 2 and taking out derivatives. Prentice-Hall, Englewood Cliffs, NJ.
k = 1 we get: |det A (1, 20)| = 3 × 10−47, |detA (1, 30)| =
4 × 10−117 , |det A (1, 40)| = 3 × 10−217 . The wave num- Chen, J. T.; Chen, I. L.; Chen, K. H.; Lee, Y. T.; Yeh,
ber k = 1 is not the eigenvalue of the problem. This is the Y. (2004): A meshless method for free vibration analysis
”background” value between extremums and one looks of circular and rectangular clamped plates using radial
for the minima of |detA (k, N)| on such background. So, basis function. Eng. Anal. Bound. Elem., vol. 28, pp.
using this technique one operates with values of the order 535–545.
∼ 10−50 − 10−500 , see [Alves and Antunes (2005); Chen,
Chen, J. T.; Chen, I. L.; Lee, Y. T. (2005): Eigensolu-
Chen, and Lee (2005)] for more detailed information.
tions of multiply connected membranes using the method
At the same time let us calculate the norm function of fundamental solutions. Eng. Anal. Bound. Elem., vol.
F(k, N) which is used to obtain the eigenvalues in the 29, pp. 166–174.
method presented. We get for ε = 0.0001: F(1, 20) =
2.13 × 10−5, F(1, 30) = 2.13 × 10−5, F(1, 40) = 2.13 × Chen, J. T.; Lin, J. H.; Kuo, S. R.; Chyuan, S. W.
10−5 . We present the values of the norm function F (k) (2001): Boundary element analysis for the Helmholtz
when k is close to eigenvalue in Tab. 14. eigenvalues problems with a multiply connected domain.
Here the number of the sources is fixed N = 30 and the Proc. R. Soc. Lond. A, vol. 457, pp. 2521–2546.
smoothing parameter ε is varied. er is the relative error in
Chen, J. T.; Liu, L. W.; Hong, H. K. (2003): Spurious
determining of the approximated eigenvalue ki and F (ki )
and true eigensolutions of Helmholtz BIEs and BEMs for
denotes the value of the norm function at this approxi-
a multiply connected problem. Proc. R. Soc. Lond. A,
mated eigenvalue. So, in the framework of the method
vol. 459, pp. 1897–1924.
presented we always deal with the values which can be
handled on PC with a single precision. Chen, W. (2005): Symmetric boundary knot method.
The method is easy to program and not expensive in the Eng. Anal. Bound. Elem., vol. 26, pp. 489–494.
CPU time. The all calculations presented in the paper
were performed using 366 MHz PC. Chen, W.; Tanaka, M. (2002): A meshless, expo-
nential convergence, integration-free, and boundary-only
RBF technique. Computers and Mathematics with Ap-
References plications, vol. 43, pp. 379–391.
Alves, C. J. S.; Antunes, P. R. S. (2005): The
Courant, R. (1943): Variational methods for the so-
Method of Fundamental Solutions applied to the calcu-
lution of problems of equilibrium and vibrations. Bull.
lation of eigenfrequencies and eigenmodes of 2D simply
Amer. Math. Soc., vol. 43, pp. 1–23.
connected shapes. CMC: Computers, Materials & Con-
tinua, vol. 3, pp. 00–00. Courant, R.; Hilbert, D. (1953): Methods of Mathe-
matical Physics. Wiley Interscience, New York.
Atluri, S. (2004): The Meshless Method (MLPG)
for Domain & Boundary Discretizations. Tech Science Fairweather, G.; Karageorghis, A. (1998): The
Press, Forsyth. method of fundamental solutions for elliptic boundary
 188 Copyright

c 2005 Tech Science Press CMC, vol.2, no.3, pp.177-188, 2005

value problems. Advances in Comp. Math., vol. 9, pp. Reutskiy, S. Y. (2002): A boundary method of Tre-
69–95. fftz type with approximate trial functions. Engineering
Analysis with Boundary Elements, vol. 26, pp. 341–353.
Golberg, M. A.; Chen, C. S. (1997): Discrete Pro-
jection Methods for Integral Equations. Computational Strang, G. (1976): Linear Algebra and Its Applica-
Mechanics Publications, Southampton. tions. Academic Press.

Golberg, M. A.; Chen, C. S. (1998): The Method Strang, G.; Fix, G. (1973): An Analysis of the Finite
of Fundamental Solutions for Potential, Helmholtz and Element Method. Prentice-Hall.
Diffusion Problems, pp. 103–176. WIT Press, 1998.
Tsangaris, T.; Smyrlis, Y. S.; Karageorghis, A. (2004):
Golub, G. H.; Loan, C. F. V. (1996): Matrix Compu- A Matrix Decomposition MFS Algorithm for Bihar-
tations. Johns Hopkins University Press, Baltimore. monic Problems in Annular Domains. CMC: Comput-
ers, Materials, & Continua, vol. 1, pp. 245–258.
Grannell, J. J.; Atluri, S. N. (1978): Boundary Meth-
ods (BEM) and Combination of BEM-FEM. GIT-ESM- Vekua, I. (1957): New Methods for Solving Elliptic
SNA-78-16. Equations. North-Holland, Amsterdam.

Hafner, C. (1990): The Generalized Multipole

Technique for Computational Electromagnetics. Artech
House Books, Boston.

Han, Z. D.; Atluri, S. N. (2003): Truly Meshless Local

Petrov-Galerkin (MLPG) Solutions of Traction & Dis-
placement BIEs. CMES: Computer Modeling in Engi-
neering & Sciences, vol. 4, pp. 665–678.

Han, Z. D.; Atluri, S. N. (2004): A Meshless Local

Petrov-Galerkin (MLPG) Approach for 3-Dimensional
Elasto-dynamics. CMC: Computers, Materials, and
Contniua, vol. 1, pp. 129–140.

Kang, S. W.; Lee, J. M. (2001): Free vibration anal-

ysis of arbitrary shaped plates with clamped edges using
wave-type functions. J. Sound. Vibr., vol. 242, pp. 9–26.

Karageorghis, A. (2001): The method of fundamen-

tal solutions for the calculation of the eigenvalues of the
Helmholtz equation. Applied Math. Letters., vol. 14, pp.
837–842.

Kupradze, V. D. (1967): On approximate solution of

problems in mathematical physics. Russian Math. Sur-
veys, vol. 22, pp. 58–108.

Morse, P. M.; Feshbach, H. (1953): Methods of The-

oretical Physics. McGraw–Hill.

Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flan-

nery, B. P. (2002): Numerical Recipescipes in C++.
Cambridge University Press.

View publication stats