- Email: [email protected]
A boundary element scheme for diffusion problems using compactly supported radial basis functions
Engineering Analysis with Boundary Elements 23 (1999) 201–209
A boundary element scheme for diffusion problems using compactly supported radial basis functions M. Zerroukat* Wessex Institute of Technology, Ashurst Lodge, Ashurst, Hampshire S040 7AA, UK Received 11 March 1998; accepted 1 October 1998
Abstract The solution of the diffusion equation using the boundary element method has a high computational cost due to the inherent time history constraint in the integral representation. This time-history dependence become impractical for problems where computations are to be performed for extended times. In general, the computation of the solution at n domain points using m boundary points at the kth time-step requires an amount of computer operations of the order O(km 2 ⫹ knm). This paper presents a scheme that requires a computational cost of the order of only O(m 2 ⫹ nm), where the dependence from the past k-steps is removed. In other words, the computational process is reduced to that in finite differences and finite elements, where the solution at every time step is dependent from that of the previous step only. The scheme uses the time dependent fundamental solution but the time integration is performed over one time-step only and the rest of the history integral is converted to a domain integral which is approximated using compactly supported radial basis functions. 䉷 1999 Elsevier Science Ltd. All rights reserved.
1. Introduction The numerical solution of time-dependent diffusion equation via the boundary integral approach have been greatly hampered by the time dependence in the boundary integral formulation. For instance, evaluating the solution require more and more computational work as time progresses. Many authors have attempted to overcome this major drawback in the solution of potential problems using the boundary element method. For instance, Demirel and Wang [1] subdivide the whole time-integral into near-history and farhistory integrals, where the near-history is evaluated in the normal convolution procedure while the far-history one is approximated, due to the fact that the fundamental solution decays with time. Davey and Hinduja [2] evaluates the nearhistory integral in the normal manner, while replacing the far-history one with a domain integral, which is approximated using weights deduced from the boundary solution at the previous timesteps. Davey and Bounds [3] used an approach whereby the fundamental solution at any time step is approximated by a linear combination of solutions evaluated at previous time steps. This procedure is a variant of that reported by Davey and Hinduja [2] and; therefore; * Tel.: ⫹ 44-1703-293223;Fax: ⫹ 44-1703-292853; e-mail: zerr@ wessex.ac.uk.
suffers from the same drawbacks. All these methods [1–3] use a truncation, where the accuracy depends on its length, which is effective only after a large number of time steps. Further, approximating a domain integral using the boundary solution only is bound to introduce errors if the solution is not smooth enough (presence of sharp and regional variations). Therefore, the truncation doesn’t really overcome the problem of time history altogether but reduces its scale slightly. Greengard and Strain [4] instead of using the free-space Green’s function, as a weighting function, they used a complementary representation obtained by Fourier series and the method of images. Although, the representation itself does not eliminate the time history, the inherent fact that the Fourier series coefficients, at any time step, can be recursively computed from those at previous steps is enough to make some savings in cpu-time. Another approach is the dual reciprocity [5], which uses a time-independent fundamental solution and treats the time-derivative as a non-homogenous term. Time stepping is carried out in similar manner as in finite differences and finite elements, and, therefore, it lacks the high accuracy associated with time-dependent fundamental solutions, especially for large time-steps [6]. Laplace or Fourier transforms can also be used to eliminate the time history constraint [7, 8]. Although very effective when the solution is sought at some instant in time without the need for the previous history, the
0955-7997/99/$ - see front matter 䉷 1999 Elsevier Science Ltd. All rights reserved. PII: S0955-799 7(98)00089-7
202
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
transformation and inversion process from real to transformed spaces is costly if the evolution of the history of the solution is sought in a time-stepping manner. Zerroukat [9] used a scheme that employs the timedependent fundamental solution in the integral representation but the time integration is performed over one time step only. The rest of the time-history integral is converted to a domain integral, which is approximated using global radial basis functions. This results in huge savings in terms of cputime and put the scheme on the same playing field as finite differences and finite elements where the solution at every time step is computed from that at the previous one only. When using global radial functions, the resulting interpolation matrix is fully populated. The system may also become ill-conditioned if very smooth radial functions are used on dense and large number of points. Although pre-conditioning methods [10] can be used to deal with this problem, both the cpu-time and ill-conditioning increase with increasing the number of collocation points. This limits the use of global radial basis functions to a maximum number of collocation points, which depends on the power of the computational platform. Other alternatives are domain decomposition [11] and compactly support radial basis functions [12]. This results in a sparse matrix coefficient and bypass the above mentioned problems related to global radial basis functions. In this paper, a scheme similar to that given previously by the same author [9], but with compactly supported radial basis functions is presented and tested herein. The scheme approximates the time-dependent domain integral as a linear combination of compactly supported radial basis functions, where the coefficients at every time step are computed from those at the previous one. This results in a sparse matrices, which can be efficiently inverted using iterative methods, hence improving even further the huge savings in cpu-time already achieved with global radial basis functions.
2. Integral approximation using radial basis functions In this section a scheme using BEM for the solution of diffusion type boundary value problems, whereby the time integration is performed over one time step only and approximating the domain integral, will be presented. The scheme uses the time-dependent fundamental solution in the integral representation and approximates the evolution of the domain integral using compactly supported radial basis functions scheme. Consider the general diffusion equation:
2u x; t k7 2 u x; t; x 僆 V 傺 R3 ; t ⬎ 0 2t
corresponding to Eq. (1) over the entire space-time domain can be written as [13]: c ju j; t k
⫹
t0
Z
dt
Z V
G
uⴱ
2u 2uⴱ ⫺ dG 2ns 2ns
u ⴱ j; x; t; t0 u x; t0 dV
2
where u* is the free space Green’s function given by: ( ) 1 r2 exp ⫺ u ⴱ j; x; t; t 4k t ⫺ t 4pk t ⫺ td=2
3
where d is the dimension of the problem, r 储j-x储 is the Euclidian distance between the field point x and the source point j, ns is the outward normal to the boundary G, and c(j) is a constant which depends on the location of the source point j and the local geometry (e.g. for a smooth boundary c(j) 1/2). Assuming that t0 t n and t t n⫹1 t n ⫹ dt, Eq. (2) is written as
c ju j; tn⫹1 k
Zt n ⫹ 1
⫹
tn
Z V
dt
Z G
uⴱ
2u 2uⴱ ⫺ dG 2ns 2ns
g j; x;dtu x; tn dV
4
where dt is the time step size and g j; x; dt a=pd=2 exp{ ⫺ a储j ⫺ x储2 }
5
where a 1/4kdt. Let p n(j,t) be defined as p n j; t
Z V
g j; x; tu x; tn dV
6
Assuming that, given a total of N boundary and domain collocation points, then u(x,t n) can be approximated by: u x; tn ⯝
N X
lnj w x; xj
7
j1
where the compactly supported radial basis function w(x,xj) is given by [12], [14]:
1
with certain conditions on the boundary G 2V, where u(x,t) denote concentration at the spatial position x at time t, 7 the gradient differential operator; V is a bounded domain in R 3, and k, the diffusivity. The integral equation
Zt
(
w x; xj
f x; xj for 0
for
0 ⱕ 储x ⫺ xj 储 ⱕ c xj 储x ⫺ xj 储 ⬎ c xj
8
where c(xj) is the compact support of w(x,xj) at xj. Using Eq.
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
be rewritten, in general, as: Z 4ij v ji ; xj ; dt c ji ; x; xj ; dtdV j
(7) and Eq. (8), then Eq. (6) can be approximated by: pn j; ddt ⯝
N X
lnj
j1
N X
lnj
Z V
Vj
j1
N X j1
N X
lnj
g j; x; dtw x; xj d V
V j
!
Z
rij
g j; x; dtf x; xj d Vj !
Z Vj
c j; x; xj ; dtdVj
lnj v j; xj ; dt {ln }T {vj }
9
j1
Z Vj
c ji ; x; xj ; dtdVj ; V j 傺 Vj Vj 傼 V
Vj
where {vj } [vj1…vjN] , vj j v j; xj ; dt R n T n n Vj c j; x; xj ; dtdVj ; {l } l1 …lN ; c j; x; xj ; dt g j; x; dtf x; xj and Vj is the support domain at xj (i.e. a distance of length 2c(xj), a circle of radius c(xj) or a sphere of radius c(xj) with xj at the centre, for one-, two- or threedimensional problems, respectively; Vj {x,0 ⱕ 储x-xj储 ⱕ c(xj)}). The notation {…} and […] are used for vector and matrix, respectively; […] T denotes matrix transpose. Applying Eq. (4) to every boundary point j j1,…,jk the following system of equations can be obtained: H{un⫹1 } ⫹ G{qn⫹1 } {pn }
10
where {un⫹1 } u j1 ; tn⫹1 …u jk ; tn⫹1 T ; {qn⫹1 } 2u=2ns j1 ; tn⫹1 …2u=2ns jk ; tn⫹1 T and {pn } pn j1 ; dt…pn jk ; dtT which is given by
411
6 6 . {pn } 6 .. 4
4k1
12
where 0 ⬍ rij ⱕ 1, a positive geometry factor given by: 8 1 when Vj 傼 V Vj > > > > > > ⬍ 1 when Vj 傼 V V j Z > > > c
ji ; x; xj ; dtdVj :
13
T
2
203
… 41N 3 7 .. 7 n {l } 4{ln } ] . 7 5 … 4kN
11
where 4ij v(ji,xj,dt). For the coefficients matrices H and G, see Wrobel and Brebbia [13] for details. It has to be noted that the entries of [4] depend on only dt and the locations of both the boundary and the internal points. As it is a common practice to keep the same collocation points and a constant time step throughout the computational process, the computation of [4] is performed only once. This reduces the computation of {p n} at every time step to simply compute the vector {l n}. The entries of [4] consist of the integration of a known function c, with variable parameters ji1,k and xi1,N, over a simple support domain Vj. For the collocation points m where the support domain Vm intersects the boundary G 2V of the domain V, i.e. part of Vm lay outside V, 4im can be calculated using techniques such as form-factors as used in the radiosity method of computer graphics to account for these partial overlaps. Denoting by V j the part of the support domain Vj that lay inside V (i.e. V j Vj 傼 V), then 4ij in Eq. (11) can
rij can also be approximated simply as rij 艑 rj 艑 V j /Vj. Since rij are geometry dependent coefficients, i.e. dependent on location of centres with respect to the domain V, they can be computed beforehand and only once. Although this will involve some extra computation, it will have little effect on the efficiency of overall time-marching scheme. In order to compute {p n} from {p n-1}, the approach reported by Kansa [16] and Zerroukat et al. [17] is adopted. Since Eq. (10) consists of computing the solution [u n⫹1 q n⫹1] T at the level time (n ⫹ 1) from those at level n, and therefore the solutions {u n} and {u n-1} at the levels n and (n1), respectively, are already known. Hence, the evaluation of {p n} from {p n-1}, i.e. computing {l n} from {l n-1}, is similar to solving a local Dirichlet problem, for t n-1 ⬍ t ⬍ t n, defined by:
2u x; t k72 x; t; x 僆 V 傺 R3 ; t 僆 tn⫺1 ; tn 2t
14
u x; t f x; t; x 僆 G 2V; t 僆 tn⫺1 ; tn
15
u x; t un⫺1 x; t tn⫺1
16
where f(x,t) and u n-1(x) are known functions. It has to be emphasised that the boundary and initial conditions Eqs. (15)-(16) are local, t n-1 ⬍ t ⬍ t n, and do not represent the boundary and initial conditions of the global diffusion problem, for which the BEM is used. The global boundary and initial conditions are treated in the same manner as in standard BEM [13]. For instance, f in Eq. (15) is a simple linear function of the boundary solutions at the time steps (n) and (n-1), e.g. f(x,t) j1(t)u n(x) ⫹ j2(t)u n-1(x). First, let us discretize Eq. (14), according to the uweighted scheme giving u x; tn ⫺ u x; tn⫺1 dtk{u72 u x; tn ⫹ 1 ⫺ u72 u x; tn⫺1 }
17
where 0 ⬍ u ⬍ 1, and t n t n-1 ⫹ dt. Using the notation u n
204
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
Table 1 Compactly supported radial basis functions for d 1,5,…,n
xj) is given by:
d1
f(r) (1-r) f(r) (1-r) 3(3r ⫹ 1) f(r) (1-r) 5(8r 2 ⫹ 5r ⫹ 1)
w xi ; xj w rij
d3
f(r) (1-r) f(r) (1-r) 4(4r ⫹ 1) f(r) (1-r) 6(35r 2 ⫹ 18r ⫹ 3)
C C2 C4
d5
f(r) (1-r) 3 f(r) (1-r) 5(5r ⫹ 1) f(r) (1-r) 7(16r 2 ⫹ 7r ⫹ 1)
C0 C2 C4
2
C0 C2 C4
lX ⫹ 2k
u(x,t ), Eq. (17) can be rearranged as: 18
where a -kudt and b kdt(1-u). Applying Eq. (7) to every collocation points, xi1,…,N, the following system is obtained: un xi ; tn u xi ; yi ; zi ; tn ⯝
0
for
rij ⬎ 1
20
alj;k r j
21
j0
n
N X
f rij for 0 ⱕ rij ⱕ 1
where rij 储xi-xj储/cj, and cj c(xj) is the compact support of w at xj (for the general case of variable compact support). f(rij) is a positive definite radial function, which according to Wendland [14], can be computed for any dimension d and a required smoothness 2k using (subscripts i, j are dropped from r for neatness):
0
f r
un ⫹ a72 un un⫺1 ⫹ b72 un⫺1
(
lnj w xi ; xj ; i 1; …; N 19
with l bd/2c ⫹ k ⫹ 1 (byc denotes the largest integer less than or equal to y) and the coefficients alj;k can be computed recursively for any 0 ⱕ s ⱕ k-1 using: alj;0 ⫺1j al0;s⫹1
where (xi, yi, zi) denotes spatial coordinates of xi, and w(xi,
lX ⫹ 2s j0
j1
alj;s⫹1
l! ; 0ⱕjⱕl l ⫺ j!j! alj;s ; al 0; s ⱖ 0 j ⫹ 2 1;s⫹1
22
1 l ⫺ ; s ⱖ 0; 2 ⱕ j ⱕ l ⫹ 2s ⫹ 2 a j j⫺2;s
Fig. 1. Comparison of CSRBF-BEM and DBEM for t 0.2, 0.3 and 0.5. N 16 (quasi-random), dt 0.1, k 0.1, m 1.0, c(x) 0.5 (constant).
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
205
Fig. 2. Typical speed-up factor SF versus dimension of sparse systems N. SF (cpu-time direct solver)/(cpu-time iterative solver). The cpu-times used to compute SF are specific to Gateway 2000 P5-120.
For convenience, let us reproduce the table of functions given by Wendland [14] for the obvious cases of d 1,3,5: It can be seen from Table 1 that only odd ds are tabulated. Every basis function that is positive definite on R d is also positive definite on R k with k ⱕ d. The entries for the space dimension d in Table 1 give the maximum possible space dimension where the basis function is positive definite. However, if for instance d 2, functions that are positive definite on R dⱖ3 can be used. The functions given for a space d are also optimal for k ⱖ d. The construction of functions using Eq. (21) leads to the same functions for a space dimension d 2n and d 2n ⫹ 1. Since the diffusion equation involves second-order derivatives, a minimum of C 2 is required for f. In general Eq. (19) can also be written with an additional polynomial on the right-hand side, see Zerroukat [9] for details. However, in this paper and for simplicity, the polynomial is omitted due to its negligible effect on the accuracy of the scheme. Rewriting Eq. (19) in matrix form, viz:
given by: 2 w11 6 6 A 6 ... 4
wN1
… .. .
w1N
3
7 .. 7 . 7 5
24
… wNN
Assuming that there are NV ⬍ N internal (domain) points and NG (N-NV) boundary points, i.e. N (NV ⫹ NG), then the (N × N) matrix A can be split into: A AV ⫹ AG, where AV wij for 1 ⱕ i ⱕ NV ; 1 ⱕ j ⱕ N; 0 elsewhere AG wij for NV ⬍ i ⱕ N; 1 ⱕ j ⱕ N; 0 elsewhere 25 Using the notation LA to denote the matrix of the same dimension as A and containing the elements w~ ij , where w~ ij Lwij, where L is a linear operator, then Eq. (18) together with Eq. (15) can be written, in a matrix form, as: {B}ln Cb {ln⫺1 } ⫹ {Fn }
26
where un Aln
23
where {u n} [un1 …unN ] T, {l n} [ln1 …lnN ] T and A is
B Ca ⫹ AG and Cr AV ⫹ r72 AV ; r a; b
27
and {F n} [0…0 fNnV ⫹1 fNnV ⫹2 …fNn ] T. Rewriting Eq. (26) and
206
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
P Fig. 3. Influence of the size of the compact support c on the average error. Eav [1/(N-2)] N⫺2 i2 兩uNu(xi,t)-uAn(xi,t)兩 (the two boundary points are not included in Eav), where uNu(xi,t) and uAn(xi,t) denote the computed and the analytical value of u at the domain point i at time t 0.5. N 12 (regular grid), dt 0.1, k 0.1, m 1.0.
[15]. This results in a further reduction in the overall cputime of the scheme.
putting it in a simpler form, {ln } D{ln⫺1 } ⫹ {En }
28 3. Numerical example
where D B⫺1 Cb and {En } B⫺1 {Fn }
29
Although B and D, in Eq. (29), look complicated, their computation is a simple and straightforward operation. Further, if the same collocation points and a constant time stepping scheme are used throughout the computational process, B and D are computed only once, hence computing {l n} from {l n-1} is a simple matrix–vector product operation of order O(N). Although Eq. (28) is valid for any u 僆 [0,1], the value of u 1/2 is used (i.e. the Crank–Nicholson scheme), hence Ca 2AV-Cb, where b kdt/2. Knowing the initial distribution u(x,t0), {l 0} can be obtained from Eq. (23), i.e. {l 0} A -1{u 0}. The computation of the solution at any time step (n ⫹ 1) involves solving Eq. (10), with the right-hand side calculated from Eq. (11), where the vector {l n} is calculated from {l n-1} using Eq. (28). Due to the sparse nature of the matrices involved, instead of using a direct solver, an iterative one based on the biconjugate gradient method (BGM) is used to solve the sparse systems
For simplicity, the present method will be referred to as CSRW-BEM (compactly supported radial basis functions boundary element method) while DBEM will refer to the direct boundary element method, which contains the original time-history integral. For validation, a simple onedimensional diffusion problem is considered. However, the conclusions reached for this simple case is expected to remain valid for other more complicated problems in higher dimensions. Let us consider the following boundary value problem governed by the following equations:
2u 22 u k 2 ; 0 ⬍ x ⬍ 1; t ⬎ 0 2t 2x
30
u 0; t u 1; t 0; t ⱖ 0
31
( u x; 0
mx;
0 ⱕ x ⱕ 1=2
m 1 ⫺ x; 1=2 ⱕ x ⱕ 1
32
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
207
√ Fig. 4. Influence of the size of the compact support (c/rmax c/ 2) on the condition number of the interpolation matrix A using 81 points in a unit square and 2 two radial basis functions f(r) (1-r) (continuous line) and f(r) (4r ⫹ 1)(1-r) 4 (dashed line).
The analytical solution is given by [18, 19]: ∞ 4m X 1 npx 2 sin exp⫺kt np u x; t 2 2 p n1 n2
33
Fig. 1 shows the computed distribution u(x,t) at different times t, using the proposed method (CSRBF-BEM), the direct method (DBEM), and the analytical solution. It can be seen that both CSRBF-BEM and DBEM are in good agreement with the analytical solution. It has already been reported in [9] that this time-marching scheme, coupled with globally supported radial basis functions, can achieve cpu-time gains of the order O(k) at every time step and O(k/ 2) for the accumulated gains, when compared to DBEM. This gain is primarily due to the fact that, unlike DBEM, the present scheme removes the time history dependence and renders the computational process similar to that in FDM and FEM, where the solution at every time step is dependent from that at the previous one only. These characteristics of the scheme in [9] are also preserved for CSRBF-BEM, but with further gains due to the use of compact instead of global supports. This results in sparse matrices instead of fully populated ones in [9], which can be solved more efficiently using iterative instead of direct solvers. Fig. 2 shows the variation of the speed-up factor SF, due the use of biconjugate-gradient instead of Gaussian
elimination, as function of the number of collocation points N. It can be seen that as N increases, the gain in cpu-time using iterative, instead of direct, solvers became substantial. The cpu-time gain is in the order of O(N 3) . This allow the possibility of using large number of points without any substantial increase in the cpu-time. Fig. 3 shows the influence of the size of the compact support on the overall accuracy of the scheme, where the accuracy increases as c increases with a concomitant increase in the ill-conditioning of the matrices. As more distant points are covered by the centre associated cloud (a synonym term of compact support widely used in meshless or point based methods), the ill-conditioning of the matrices involved increases. The ill-conditioning increases not only with the size of the compact support, but also with the smoothness of the interpolant (radial basis functions). This is illustrated in Fig. 4, where typical condition numbers for the interpolation matrix A are computed using 81 points/centres in a unit square [0,1] 2, and two radial functions f(r) (1-r) 2 and f(r) (4r ⫹ 1)(1-r) 4 with C 0 and C 2 continuity, respectively. It can be seen that the condition number of the interpolation matrix A increases rapidly with the size of the compact support and the smoothness of the radial function used. Further, tests were carried out with variable, adaptive compact support c(x), and the results suggest that the
208
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
Table 2 Comparison of solutions due to DBEM and CSRBF-BEM with variable compact support c(x). u(x,t) at different positions x at different times t. k 0.05, m 0.5, dt 0.1, N 15 (regular grid). Note: For every time t, the first, second, and third row correspond to u(x,t) computed using DBEM, CSRBF-BEM, and the analytical solution, respectively. The compact support was taken as: c(x) 0.4x ⫹ 0.45 for 0 ⱕ x ⱕ 0.5, and c(x) -0.4x ⫹ 0.85 for 0.5 ⱕ x ⱕ 1, i.e. 0.45 ⱕ c(x) ⱕ 0.65 Methods
x 0.07143
0.21429
0.35714
0.50000
0.64286
t 0.2 DBEM CSRBF ANAL
0.03573 0.02889 0.03567
0.10585 0.10133 0.10601
0.16605 0.16731 0.16701
0.19197 0.19346 0.19358
0.16601 0.16659 0.16701
t 0.4 DBEM CSRBF ANAL
0.03466 0.03107 0.03469
0.09990 0.09691 0.10027
0.14981 0.14993 0.15068
0.16906 0.16903 0.17021
0.14973 0.14857 0.15068
t 0.6 DBEM CSRBF ANAL
0.03251 0.03020 0.03256
0.09214 0.08980 0.09255
0.13536 0.13455 0.13613
0.15133 0.15026 0.15228
0.13529 0.13344 0.13613
t 0.8 DBEM CSRBF ANAL
0.02994 0.02810 0.02998
0.08416 0.08187 0.08455
0.12248 0.12095 0.12317
0.13637 0.13458 0.13719
0.12242 0.12011 0.12317
t 1.0 DBEM CSRBF ANAL
methods when it comes to computational cost. The local nature of compactly supported radial basis functions, gives the scheme some flexibility for heterogeneous mediums. The resulting sparse matrices also allow the problem to be solved using a sufficiently large number of points without a substantial penalty on cpu-time. Adaptive points density and compact support can also be easily incorporated to deal with regional variations in the smoothness of the solution and geometrical complexities. Further, unlike element based integration, the radial-basis function integration scheme is remarkably simple, especially in higher dimensions.
Acknowledgements This research project is supported by the Department of Trade and Industry of the United Kingdom, and forms part of the action COST-512 of the European Commission. I also wish to thank Prof. R. Schaback and Dr Holger Wendland of the Institute for Numerical and Applied Mathematics, University of Goettingen, Germany, for their valuable discussions on many occasions.
References 0.02733 0.02565 0.02736
0.07653 0.07410 0.07689
0.11091 0.10882 0.11152
0.12325 0.12093 0.12398
0.11087 0.10819 0.11152
accuracy can be tuned accordingly, by varying the compact support. Compact support can be increased in regions of sharp variations or increased geometrical complexity, whereas decreased elsewhere. In Table 2, a comparison of DBEM, analytical, and CSRBF-BEM with a variable c(x) is shown. For the example in Table 2, a simple linear relation is used for c(x). Due to the fact that the number of points falling within a constant support domain decreases as we move towards the boundaries, c(x) is chosen in a way to compensate for the latter by increasing c(x) as we move towards the boundaries. The numerical results are in a good agreement with DBEM and analytical solutions.
4. Conclusions The numerical results clearly demonstrate that despite the fact that the direct boundary element method is in general more accurate, the present scheme achieves a relatively comparable accuracy but with huge savings in the computational cost. Unlike the direct method, the present scheme removes the time-history dependence and allows the solution at any time step to be computed with a relatively constant cpu-time. This put the scheme on the same level playing field with finite-differences and finite elements
[1] Demirel V, Wang S. An efficient boundary element method for twodimensional transient wave propagation problems. Appl Math Modelling 1987;11:411–416. [2] Davey K, Hinduja S. An improved procedure for solving transient heat conduction problems using the boundary element method. Int J Numer Meth Engng 1989;28:2293–2306. [3] Davey K, Bounds S. Source-weighted domain integral approximation for linear transient heat conduction. Int J Numer Meth Engng 1996;39:1775–1790. [4] Greengard L, Strain J. A fast algorithm for the evaluation of heat potentials. Comm Pure Appl Math 1990;43:949–963. [5] Partridge PA, Brebbia CA, Wrobel LC. Dual reciprocity boundary element method. Southampton and London: Computational Mechanics Publications and Elsevier, 1992. [6] Zerroukat M, Wrobel LC. A boundary element method for multiple moving boundary problems. J Comput Phys 1997;138:501–519. [7] Moridis GJ, Reddell DL. The Laplace transform boundary element method for the solution of diffusion type equations. In: Brebbia CA, Gipson GS, editors. Boundary elements XIII. Southampton and London: Computational Mechanics Publications and Elsevier, 1991. pp. 83-97. [8] Hill LR, Faxris TN. Fast Fourier transform of spectral boundary elements for transient heat conduction. Int J Numer Meth Heat Fluid Flow 1995;57:813–827. [9] Zerroukat M. A fast boundary element algorithm for time-dependent potential problems. Appl Math Modelling 1998;22:183–196. [10] Dyn N, Levin D, Rippa S. Numerical procedures for global surface fitting of scattered data by radial functions. SIAM J Sci Stat Comput 1986;7:639–659. [11] Dubal MR. Domain decomposition and local refinement for multiquadric approximations. J Appl Sci Comput 1994;1:146–171. [12] Floater MS, Iske A. Multistep scattered data interpolation using compactly supported radial basis functions. J Comput Appl Math 1996;73:65–78. [13] Wrobel LC, Brebbia CA. Boundary element methods in heat transfer.
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209 Southampton and London: Computational Mechanics Publications and Elsevier, 1992. [14] Wendland H. Error estimates for interpolation by compactly supported radial basis functions of minimal degree. J Approx Theory 1998;93:258-272. (Reprints in: http://www.num.math.uni-goettingen.de/wendland/papers.html) [15] Press WH, Teukolsky SA, Vetterling WT, Flannery BR. Numerical recipes in Fortran, 2nd ed. Cambridge: Cambridge University Press, 1992. [16] Kansa EJ. Multiquadrics: a scattered data approximation scheme with
209
applications to computational fluid-dynamics. Comput Math Appl 1990;19:147–161. [17] Zerroukat M, Power H, Chen CS. A numerical method for heat transfer problems using collocation and radial basis functions. Int J Numer Meth Engng 1998;42:1263–1278. [18] Carslaw HS, Jaeger JC. Conduction of heat in solids, 2nd ed. Oxford: Clarendon Press, 1959. [19] Haberland C, Lahrmann A. A comparative investigation on recurrence formulae in finite differences methods. Int J Numer Meth Engng 1988;25:593–609.
A boundary element scheme for diffusion problems using compactly supported radial basis functions M. Zerroukat* Wessex Institute of Technology, Ashurst Lodge, Ashurst, Hampshire S040 7AA, UK Received 11 March 1998; accepted 1 October 1998
Abstract The solution of the diffusion equation using the boundary element method has a high computational cost due to the inherent time history constraint in the integral representation. This time-history dependence become impractical for problems where computations are to be performed for extended times. In general, the computation of the solution at n domain points using m boundary points at the kth time-step requires an amount of computer operations of the order O(km 2 ⫹ knm). This paper presents a scheme that requires a computational cost of the order of only O(m 2 ⫹ nm), where the dependence from the past k-steps is removed. In other words, the computational process is reduced to that in finite differences and finite elements, where the solution at every time step is dependent from that of the previous step only. The scheme uses the time dependent fundamental solution but the time integration is performed over one time-step only and the rest of the history integral is converted to a domain integral which is approximated using compactly supported radial basis functions. 䉷 1999 Elsevier Science Ltd. All rights reserved.
1. Introduction The numerical solution of time-dependent diffusion equation via the boundary integral approach have been greatly hampered by the time dependence in the boundary integral formulation. For instance, evaluating the solution require more and more computational work as time progresses. Many authors have attempted to overcome this major drawback in the solution of potential problems using the boundary element method. For instance, Demirel and Wang [1] subdivide the whole time-integral into near-history and farhistory integrals, where the near-history is evaluated in the normal convolution procedure while the far-history one is approximated, due to the fact that the fundamental solution decays with time. Davey and Hinduja [2] evaluates the nearhistory integral in the normal manner, while replacing the far-history one with a domain integral, which is approximated using weights deduced from the boundary solution at the previous timesteps. Davey and Bounds [3] used an approach whereby the fundamental solution at any time step is approximated by a linear combination of solutions evaluated at previous time steps. This procedure is a variant of that reported by Davey and Hinduja [2] and; therefore; * Tel.: ⫹ 44-1703-293223;Fax: ⫹ 44-1703-292853; e-mail: zerr@ wessex.ac.uk.
suffers from the same drawbacks. All these methods [1–3] use a truncation, where the accuracy depends on its length, which is effective only after a large number of time steps. Further, approximating a domain integral using the boundary solution only is bound to introduce errors if the solution is not smooth enough (presence of sharp and regional variations). Therefore, the truncation doesn’t really overcome the problem of time history altogether but reduces its scale slightly. Greengard and Strain [4] instead of using the free-space Green’s function, as a weighting function, they used a complementary representation obtained by Fourier series and the method of images. Although, the representation itself does not eliminate the time history, the inherent fact that the Fourier series coefficients, at any time step, can be recursively computed from those at previous steps is enough to make some savings in cpu-time. Another approach is the dual reciprocity [5], which uses a time-independent fundamental solution and treats the time-derivative as a non-homogenous term. Time stepping is carried out in similar manner as in finite differences and finite elements, and, therefore, it lacks the high accuracy associated with time-dependent fundamental solutions, especially for large time-steps [6]. Laplace or Fourier transforms can also be used to eliminate the time history constraint [7, 8]. Although very effective when the solution is sought at some instant in time without the need for the previous history, the
0955-7997/99/$ - see front matter 䉷 1999 Elsevier Science Ltd. All rights reserved. PII: S0955-799 7(98)00089-7
202
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
transformation and inversion process from real to transformed spaces is costly if the evolution of the history of the solution is sought in a time-stepping manner. Zerroukat [9] used a scheme that employs the timedependent fundamental solution in the integral representation but the time integration is performed over one time step only. The rest of the time-history integral is converted to a domain integral, which is approximated using global radial basis functions. This results in huge savings in terms of cputime and put the scheme on the same playing field as finite differences and finite elements where the solution at every time step is computed from that at the previous one only. When using global radial functions, the resulting interpolation matrix is fully populated. The system may also become ill-conditioned if very smooth radial functions are used on dense and large number of points. Although pre-conditioning methods [10] can be used to deal with this problem, both the cpu-time and ill-conditioning increase with increasing the number of collocation points. This limits the use of global radial basis functions to a maximum number of collocation points, which depends on the power of the computational platform. Other alternatives are domain decomposition [11] and compactly support radial basis functions [12]. This results in a sparse matrix coefficient and bypass the above mentioned problems related to global radial basis functions. In this paper, a scheme similar to that given previously by the same author [9], but with compactly supported radial basis functions is presented and tested herein. The scheme approximates the time-dependent domain integral as a linear combination of compactly supported radial basis functions, where the coefficients at every time step are computed from those at the previous one. This results in a sparse matrices, which can be efficiently inverted using iterative methods, hence improving even further the huge savings in cpu-time already achieved with global radial basis functions.
2. Integral approximation using radial basis functions In this section a scheme using BEM for the solution of diffusion type boundary value problems, whereby the time integration is performed over one time step only and approximating the domain integral, will be presented. The scheme uses the time-dependent fundamental solution in the integral representation and approximates the evolution of the domain integral using compactly supported radial basis functions scheme. Consider the general diffusion equation:
2u x; t k7 2 u x; t; x 僆 V 傺 R3 ; t ⬎ 0 2t
corresponding to Eq. (1) over the entire space-time domain can be written as [13]: c ju j; t k
⫹
t0
Z
dt
Z V
G
uⴱ
2u 2uⴱ ⫺ dG 2ns 2ns
u ⴱ j; x; t; t0 u x; t0 dV
2
where u* is the free space Green’s function given by: ( ) 1 r2 exp ⫺ u ⴱ j; x; t; t 4k t ⫺ t 4pk t ⫺ td=2
3
where d is the dimension of the problem, r 储j-x储 is the Euclidian distance between the field point x and the source point j, ns is the outward normal to the boundary G, and c(j) is a constant which depends on the location of the source point j and the local geometry (e.g. for a smooth boundary c(j) 1/2). Assuming that t0 t n and t t n⫹1 t n ⫹ dt, Eq. (2) is written as
c ju j; tn⫹1 k
Zt n ⫹ 1
⫹
tn
Z V
dt
Z G
uⴱ
2u 2uⴱ ⫺ dG 2ns 2ns
g j; x;dtu x; tn dV
4
where dt is the time step size and g j; x; dt a=pd=2 exp{ ⫺ a储j ⫺ x储2 }
5
where a 1/4kdt. Let p n(j,t) be defined as p n j; t
Z V
g j; x; tu x; tn dV
6
Assuming that, given a total of N boundary and domain collocation points, then u(x,t n) can be approximated by: u x; tn ⯝
N X
lnj w x; xj
7
j1
where the compactly supported radial basis function w(x,xj) is given by [12], [14]:
1
with certain conditions on the boundary G 2V, where u(x,t) denote concentration at the spatial position x at time t, 7 the gradient differential operator; V is a bounded domain in R 3, and k, the diffusivity. The integral equation
Zt
(
w x; xj
f x; xj for 0
for
0 ⱕ 储x ⫺ xj 储 ⱕ c xj 储x ⫺ xj 储 ⬎ c xj
8
where c(xj) is the compact support of w(x,xj) at xj. Using Eq.
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
be rewritten, in general, as: Z 4ij v ji ; xj ; dt c ji ; x; xj ; dtdV j
(7) and Eq. (8), then Eq. (6) can be approximated by: pn j; ddt ⯝
N X
lnj
j1
N X
lnj
Z V
Vj
j1
N X j1
N X
lnj
g j; x; dtw x; xj d V
V j
!
Z
rij
g j; x; dtf x; xj d Vj !
Z Vj
c j; x; xj ; dtdVj
lnj v j; xj ; dt {ln }T {vj }
9
j1
Z Vj
c ji ; x; xj ; dtdVj ; V j 傺 Vj Vj 傼 V
Vj
where {vj } [vj1…vjN] , vj j v j; xj ; dt R n T n n Vj c j; x; xj ; dtdVj ; {l } l1 …lN ; c j; x; xj ; dt g j; x; dtf x; xj and Vj is the support domain at xj (i.e. a distance of length 2c(xj), a circle of radius c(xj) or a sphere of radius c(xj) with xj at the centre, for one-, two- or threedimensional problems, respectively; Vj {x,0 ⱕ 储x-xj储 ⱕ c(xj)}). The notation {…} and […] are used for vector and matrix, respectively; […] T denotes matrix transpose. Applying Eq. (4) to every boundary point j j1,…,jk the following system of equations can be obtained: H{un⫹1 } ⫹ G{qn⫹1 } {pn }
10
where {un⫹1 } u j1 ; tn⫹1 …u jk ; tn⫹1 T ; {qn⫹1 } 2u=2ns j1 ; tn⫹1 …2u=2ns jk ; tn⫹1 T and {pn } pn j1 ; dt…pn jk ; dtT which is given by
411
6 6 . {pn } 6 .. 4
4k1
12
where 0 ⬍ rij ⱕ 1, a positive geometry factor given by: 8 1 when Vj 傼 V Vj > > > > >
T
2
203
… 41N 3 7 .. 7 n {l } 4{ln } ] . 7 5 … 4kN
11
where 4ij v(ji,xj,dt). For the coefficients matrices H and G, see Wrobel and Brebbia [13] for details. It has to be noted that the entries of [4] depend on only dt and the locations of both the boundary and the internal points. As it is a common practice to keep the same collocation points and a constant time step throughout the computational process, the computation of [4] is performed only once. This reduces the computation of {p n} at every time step to simply compute the vector {l n}. The entries of [4] consist of the integration of a known function c, with variable parameters ji1,k and xi1,N, over a simple support domain Vj. For the collocation points m where the support domain Vm intersects the boundary G 2V of the domain V, i.e. part of Vm lay outside V, 4im can be calculated using techniques such as form-factors as used in the radiosity method of computer graphics to account for these partial overlaps. Denoting by V j the part of the support domain Vj that lay inside V (i.e. V j Vj 傼 V), then 4ij in Eq. (11) can
rij can also be approximated simply as rij 艑 rj 艑 V j /Vj. Since rij are geometry dependent coefficients, i.e. dependent on location of centres with respect to the domain V, they can be computed beforehand and only once. Although this will involve some extra computation, it will have little effect on the efficiency of overall time-marching scheme. In order to compute {p n} from {p n-1}, the approach reported by Kansa [16] and Zerroukat et al. [17] is adopted. Since Eq. (10) consists of computing the solution [u n⫹1 q n⫹1] T at the level time (n ⫹ 1) from those at level n, and therefore the solutions {u n} and {u n-1} at the levels n and (n1), respectively, are already known. Hence, the evaluation of {p n} from {p n-1}, i.e. computing {l n} from {l n-1}, is similar to solving a local Dirichlet problem, for t n-1 ⬍ t ⬍ t n, defined by:
2u x; t k72 x; t; x 僆 V 傺 R3 ; t 僆 tn⫺1 ; tn 2t
14
u x; t f x; t; x 僆 G 2V; t 僆 tn⫺1 ; tn
15
u x; t un⫺1 x; t tn⫺1
16
where f(x,t) and u n-1(x) are known functions. It has to be emphasised that the boundary and initial conditions Eqs. (15)-(16) are local, t n-1 ⬍ t ⬍ t n, and do not represent the boundary and initial conditions of the global diffusion problem, for which the BEM is used. The global boundary and initial conditions are treated in the same manner as in standard BEM [13]. For instance, f in Eq. (15) is a simple linear function of the boundary solutions at the time steps (n) and (n-1), e.g. f(x,t) j1(t)u n(x) ⫹ j2(t)u n-1(x). First, let us discretize Eq. (14), according to the uweighted scheme giving u x; tn ⫺ u x; tn⫺1 dtk{u72 u x; tn ⫹ 1 ⫺ u72 u x; tn⫺1 }
17
where 0 ⬍ u ⬍ 1, and t n t n-1 ⫹ dt. Using the notation u n
204
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
Table 1 Compactly supported radial basis functions for d 1,5,…,n
xj) is given by:
d1
f(r) (1-r) f(r) (1-r) 3(3r ⫹ 1) f(r) (1-r) 5(8r 2 ⫹ 5r ⫹ 1)
w xi ; xj w rij
d3
f(r) (1-r) f(r) (1-r) 4(4r ⫹ 1) f(r) (1-r) 6(35r 2 ⫹ 18r ⫹ 3)
C C2 C4
d5
f(r) (1-r) 3 f(r) (1-r) 5(5r ⫹ 1) f(r) (1-r) 7(16r 2 ⫹ 7r ⫹ 1)
C0 C2 C4
2
C0 C2 C4
lX ⫹ 2k
u(x,t ), Eq. (17) can be rearranged as: 18
where a -kudt and b kdt(1-u). Applying Eq. (7) to every collocation points, xi1,…,N, the following system is obtained: un xi ; tn u xi ; yi ; zi ; tn ⯝
0
for
rij ⬎ 1
20
alj;k r j
21
j0
n
N X
f rij for 0 ⱕ rij ⱕ 1
where rij 储xi-xj储/cj, and cj c(xj) is the compact support of w at xj (for the general case of variable compact support). f(rij) is a positive definite radial function, which according to Wendland [14], can be computed for any dimension d and a required smoothness 2k using (subscripts i, j are dropped from r for neatness):
0
f r
un ⫹ a72 un un⫺1 ⫹ b72 un⫺1
(
lnj w xi ; xj ; i 1; …; N 19
with l bd/2c ⫹ k ⫹ 1 (byc denotes the largest integer less than or equal to y) and the coefficients alj;k can be computed recursively for any 0 ⱕ s ⱕ k-1 using: alj;0 ⫺1j al0;s⫹1
where (xi, yi, zi) denotes spatial coordinates of xi, and w(xi,
lX ⫹ 2s j0
j1
alj;s⫹1
l! ; 0ⱕjⱕl l ⫺ j!j! alj;s ; al 0; s ⱖ 0 j ⫹ 2 1;s⫹1
22
1 l ⫺ ; s ⱖ 0; 2 ⱕ j ⱕ l ⫹ 2s ⫹ 2 a j j⫺2;s
Fig. 1. Comparison of CSRBF-BEM and DBEM for t 0.2, 0.3 and 0.5. N 16 (quasi-random), dt 0.1, k 0.1, m 1.0, c(x) 0.5 (constant).
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
205
Fig. 2. Typical speed-up factor SF versus dimension of sparse systems N. SF (cpu-time direct solver)/(cpu-time iterative solver). The cpu-times used to compute SF are specific to Gateway 2000 P5-120.
For convenience, let us reproduce the table of functions given by Wendland [14] for the obvious cases of d 1,3,5: It can be seen from Table 1 that only odd ds are tabulated. Every basis function that is positive definite on R d is also positive definite on R k with k ⱕ d. The entries for the space dimension d in Table 1 give the maximum possible space dimension where the basis function is positive definite. However, if for instance d 2, functions that are positive definite on R dⱖ3 can be used. The functions given for a space d are also optimal for k ⱖ d. The construction of functions using Eq. (21) leads to the same functions for a space dimension d 2n and d 2n ⫹ 1. Since the diffusion equation involves second-order derivatives, a minimum of C 2 is required for f. In general Eq. (19) can also be written with an additional polynomial on the right-hand side, see Zerroukat [9] for details. However, in this paper and for simplicity, the polynomial is omitted due to its negligible effect on the accuracy of the scheme. Rewriting Eq. (19) in matrix form, viz:
given by: 2 w11 6 6 A 6 ... 4
wN1
… .. .
w1N
3
7 .. 7 . 7 5
24
… wNN
Assuming that there are NV ⬍ N internal (domain) points and NG (N-NV) boundary points, i.e. N (NV ⫹ NG), then the (N × N) matrix A can be split into: A AV ⫹ AG, where AV wij for 1 ⱕ i ⱕ NV ; 1 ⱕ j ⱕ N; 0 elsewhere AG wij for NV ⬍ i ⱕ N; 1 ⱕ j ⱕ N; 0 elsewhere 25 Using the notation LA to denote the matrix of the same dimension as A and containing the elements w~ ij , where w~ ij Lwij, where L is a linear operator, then Eq. (18) together with Eq. (15) can be written, in a matrix form, as: {B}ln Cb {ln⫺1 } ⫹ {Fn }
26
where un Aln
23
where {u n} [un1 …unN ] T, {l n} [ln1 …lnN ] T and A is
B Ca ⫹ AG and Cr AV ⫹ r72 AV ; r a; b
27
and {F n} [0…0 fNnV ⫹1 fNnV ⫹2 …fNn ] T. Rewriting Eq. (26) and
206
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
P Fig. 3. Influence of the size of the compact support c on the average error. Eav [1/(N-2)] N⫺2 i2 兩uNu(xi,t)-uAn(xi,t)兩 (the two boundary points are not included in Eav), where uNu(xi,t) and uAn(xi,t) denote the computed and the analytical value of u at the domain point i at time t 0.5. N 12 (regular grid), dt 0.1, k 0.1, m 1.0.
[15]. This results in a further reduction in the overall cputime of the scheme.
putting it in a simpler form, {ln } D{ln⫺1 } ⫹ {En }
28 3. Numerical example
where D B⫺1 Cb and {En } B⫺1 {Fn }
29
Although B and D, in Eq. (29), look complicated, their computation is a simple and straightforward operation. Further, if the same collocation points and a constant time stepping scheme are used throughout the computational process, B and D are computed only once, hence computing {l n} from {l n-1} is a simple matrix–vector product operation of order O(N). Although Eq. (28) is valid for any u 僆 [0,1], the value of u 1/2 is used (i.e. the Crank–Nicholson scheme), hence Ca 2AV-Cb, where b kdt/2. Knowing the initial distribution u(x,t0), {l 0} can be obtained from Eq. (23), i.e. {l 0} A -1{u 0}. The computation of the solution at any time step (n ⫹ 1) involves solving Eq. (10), with the right-hand side calculated from Eq. (11), where the vector {l n} is calculated from {l n-1} using Eq. (28). Due to the sparse nature of the matrices involved, instead of using a direct solver, an iterative one based on the biconjugate gradient method (BGM) is used to solve the sparse systems
For simplicity, the present method will be referred to as CSRW-BEM (compactly supported radial basis functions boundary element method) while DBEM will refer to the direct boundary element method, which contains the original time-history integral. For validation, a simple onedimensional diffusion problem is considered. However, the conclusions reached for this simple case is expected to remain valid for other more complicated problems in higher dimensions. Let us consider the following boundary value problem governed by the following equations:
2u 22 u k 2 ; 0 ⬍ x ⬍ 1; t ⬎ 0 2t 2x
30
u 0; t u 1; t 0; t ⱖ 0
31
( u x; 0
mx;
0 ⱕ x ⱕ 1=2
m 1 ⫺ x; 1=2 ⱕ x ⱕ 1
32
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
207
√ Fig. 4. Influence of the size of the compact support (c/rmax c/ 2) on the condition number of the interpolation matrix A using 81 points in a unit square and 2 two radial basis functions f(r) (1-r) (continuous line) and f(r) (4r ⫹ 1)(1-r) 4 (dashed line).
The analytical solution is given by [18, 19]: ∞ 4m X 1 npx 2 sin exp⫺kt np u x; t 2 2 p n1 n2
33
Fig. 1 shows the computed distribution u(x,t) at different times t, using the proposed method (CSRBF-BEM), the direct method (DBEM), and the analytical solution. It can be seen that both CSRBF-BEM and DBEM are in good agreement with the analytical solution. It has already been reported in [9] that this time-marching scheme, coupled with globally supported radial basis functions, can achieve cpu-time gains of the order O(k) at every time step and O(k/ 2) for the accumulated gains, when compared to DBEM. This gain is primarily due to the fact that, unlike DBEM, the present scheme removes the time history dependence and renders the computational process similar to that in FDM and FEM, where the solution at every time step is dependent from that at the previous one only. These characteristics of the scheme in [9] are also preserved for CSRBF-BEM, but with further gains due to the use of compact instead of global supports. This results in sparse matrices instead of fully populated ones in [9], which can be solved more efficiently using iterative instead of direct solvers. Fig. 2 shows the variation of the speed-up factor SF, due the use of biconjugate-gradient instead of Gaussian
elimination, as function of the number of collocation points N. It can be seen that as N increases, the gain in cpu-time using iterative, instead of direct, solvers became substantial. The cpu-time gain is in the order of O(N 3) . This allow the possibility of using large number of points without any substantial increase in the cpu-time. Fig. 3 shows the influence of the size of the compact support on the overall accuracy of the scheme, where the accuracy increases as c increases with a concomitant increase in the ill-conditioning of the matrices. As more distant points are covered by the centre associated cloud (a synonym term of compact support widely used in meshless or point based methods), the ill-conditioning of the matrices involved increases. The ill-conditioning increases not only with the size of the compact support, but also with the smoothness of the interpolant (radial basis functions). This is illustrated in Fig. 4, where typical condition numbers for the interpolation matrix A are computed using 81 points/centres in a unit square [0,1] 2, and two radial functions f(r) (1-r) 2 and f(r) (4r ⫹ 1)(1-r) 4 with C 0 and C 2 continuity, respectively. It can be seen that the condition number of the interpolation matrix A increases rapidly with the size of the compact support and the smoothness of the radial function used. Further, tests were carried out with variable, adaptive compact support c(x), and the results suggest that the
208
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209
Table 2 Comparison of solutions due to DBEM and CSRBF-BEM with variable compact support c(x). u(x,t) at different positions x at different times t. k 0.05, m 0.5, dt 0.1, N 15 (regular grid). Note: For every time t, the first, second, and third row correspond to u(x,t) computed using DBEM, CSRBF-BEM, and the analytical solution, respectively. The compact support was taken as: c(x) 0.4x ⫹ 0.45 for 0 ⱕ x ⱕ 0.5, and c(x) -0.4x ⫹ 0.85 for 0.5 ⱕ x ⱕ 1, i.e. 0.45 ⱕ c(x) ⱕ 0.65 Methods
x 0.07143
0.21429
0.35714
0.50000
0.64286
t 0.2 DBEM CSRBF ANAL
0.03573 0.02889 0.03567
0.10585 0.10133 0.10601
0.16605 0.16731 0.16701
0.19197 0.19346 0.19358
0.16601 0.16659 0.16701
t 0.4 DBEM CSRBF ANAL
0.03466 0.03107 0.03469
0.09990 0.09691 0.10027
0.14981 0.14993 0.15068
0.16906 0.16903 0.17021
0.14973 0.14857 0.15068
t 0.6 DBEM CSRBF ANAL
0.03251 0.03020 0.03256
0.09214 0.08980 0.09255
0.13536 0.13455 0.13613
0.15133 0.15026 0.15228
0.13529 0.13344 0.13613
t 0.8 DBEM CSRBF ANAL
0.02994 0.02810 0.02998
0.08416 0.08187 0.08455
0.12248 0.12095 0.12317
0.13637 0.13458 0.13719
0.12242 0.12011 0.12317
t 1.0 DBEM CSRBF ANAL
methods when it comes to computational cost. The local nature of compactly supported radial basis functions, gives the scheme some flexibility for heterogeneous mediums. The resulting sparse matrices also allow the problem to be solved using a sufficiently large number of points without a substantial penalty on cpu-time. Adaptive points density and compact support can also be easily incorporated to deal with regional variations in the smoothness of the solution and geometrical complexities. Further, unlike element based integration, the radial-basis function integration scheme is remarkably simple, especially in higher dimensions.
Acknowledgements This research project is supported by the Department of Trade and Industry of the United Kingdom, and forms part of the action COST-512 of the European Commission. I also wish to thank Prof. R. Schaback and Dr Holger Wendland of the Institute for Numerical and Applied Mathematics, University of Goettingen, Germany, for their valuable discussions on many occasions.
References 0.02733 0.02565 0.02736
0.07653 0.07410 0.07689
0.11091 0.10882 0.11152
0.12325 0.12093 0.12398
0.11087 0.10819 0.11152
accuracy can be tuned accordingly, by varying the compact support. Compact support can be increased in regions of sharp variations or increased geometrical complexity, whereas decreased elsewhere. In Table 2, a comparison of DBEM, analytical, and CSRBF-BEM with a variable c(x) is shown. For the example in Table 2, a simple linear relation is used for c(x). Due to the fact that the number of points falling within a constant support domain decreases as we move towards the boundaries, c(x) is chosen in a way to compensate for the latter by increasing c(x) as we move towards the boundaries. The numerical results are in a good agreement with DBEM and analytical solutions.
4. Conclusions The numerical results clearly demonstrate that despite the fact that the direct boundary element method is in general more accurate, the present scheme achieves a relatively comparable accuracy but with huge savings in the computational cost. Unlike the direct method, the present scheme removes the time-history dependence and allows the solution at any time step to be computed with a relatively constant cpu-time. This put the scheme on the same level playing field with finite-differences and finite elements
[1] Demirel V, Wang S. An efficient boundary element method for twodimensional transient wave propagation problems. Appl Math Modelling 1987;11:411–416. [2] Davey K, Hinduja S. An improved procedure for solving transient heat conduction problems using the boundary element method. Int J Numer Meth Engng 1989;28:2293–2306. [3] Davey K, Bounds S. Source-weighted domain integral approximation for linear transient heat conduction. Int J Numer Meth Engng 1996;39:1775–1790. [4] Greengard L, Strain J. A fast algorithm for the evaluation of heat potentials. Comm Pure Appl Math 1990;43:949–963. [5] Partridge PA, Brebbia CA, Wrobel LC. Dual reciprocity boundary element method. Southampton and London: Computational Mechanics Publications and Elsevier, 1992. [6] Zerroukat M, Wrobel LC. A boundary element method for multiple moving boundary problems. J Comput Phys 1997;138:501–519. [7] Moridis GJ, Reddell DL. The Laplace transform boundary element method for the solution of diffusion type equations. In: Brebbia CA, Gipson GS, editors. Boundary elements XIII. Southampton and London: Computational Mechanics Publications and Elsevier, 1991. pp. 83-97. [8] Hill LR, Faxris TN. Fast Fourier transform of spectral boundary elements for transient heat conduction. Int J Numer Meth Heat Fluid Flow 1995;57:813–827. [9] Zerroukat M. A fast boundary element algorithm for time-dependent potential problems. Appl Math Modelling 1998;22:183–196. [10] Dyn N, Levin D, Rippa S. Numerical procedures for global surface fitting of scattered data by radial functions. SIAM J Sci Stat Comput 1986;7:639–659. [11] Dubal MR. Domain decomposition and local refinement for multiquadric approximations. J Appl Sci Comput 1994;1:146–171. [12] Floater MS, Iske A. Multistep scattered data interpolation using compactly supported radial basis functions. J Comput Appl Math 1996;73:65–78. [13] Wrobel LC, Brebbia CA. Boundary element methods in heat transfer.
M. Zerroukat / Engineering Analysis with Boundary Elements 23 (1999) 201–209 Southampton and London: Computational Mechanics Publications and Elsevier, 1992. [14] Wendland H. Error estimates for interpolation by compactly supported radial basis functions of minimal degree. J Approx Theory 1998;93:258-272. (Reprints in: http://www.num.math.uni-goettingen.de/wendland/papers.html) [15] Press WH, Teukolsky SA, Vetterling WT, Flannery BR. Numerical recipes in Fortran, 2nd ed. Cambridge: Cambridge University Press, 1992. [16] Kansa EJ. Multiquadrics: a scattered data approximation scheme with
209
applications to computational fluid-dynamics. Comput Math Appl 1990;19:147–161. [17] Zerroukat M, Power H, Chen CS. A numerical method for heat transfer problems using collocation and radial basis functions. Int J Numer Meth Engng 1998;42:1263–1278. [18] Carslaw HS, Jaeger JC. Conduction of heat in solids, 2nd ed. Oxford: Clarendon Press, 1959. [19] Haberland C, Lahrmann A. A comparative investigation on recurrence formulae in finite differences methods. Int J Numer Meth Engng 1988;25:593–609.