EngineeringAnalysis withBoundaryElements13 (1994) 387-396 1994ElsevierScienoeLimited Printed in Great Britain. All tights reserved 0955-7997/94/$07.00
ELSEVIER
Technical Note On the choice of interpolation functions used in the dual-reciprocity boundary-element method Yinglong Zhang & Songping Zhu Department of Mathematics, University of WoUongong, PO Box 1144, Wollongong, New South Wales 2500, Australia
(Received 13 October 1993; revised version received 13 April 1994; accepted 26 April 1994) The dual-reciprocity boundary-element method is a very powerful technique for solving general elliptic equations of the type V2u -- b. In this method, a series of interpolation functions is used to approximate b in order to convert the associated domain integral, which it is necessary to evaluate in a traditional boundaryelement analysis,into boundary integralsonly. Hence the choice of interpolation functions has direct effectson the numerical results.According to Partridge and Brebbia, the adoption of a comparatively simple form of interpolationfunction gives the best results.Unfortunately, when b contains partial derivativesof the unknown function u(x,y), the adoption of such a type of interpolationfunction inevitablyleads to the creation of singularitieson allboundary and internalnodes used in a dual-reciprocityboundary-element analysis,as was pointed out by Zhu and Zhang in 1992. To avoid this problem, a functional transformation, which applies only to linear governing equations, can be employed to eliminate these derivative terms and thus to obtain better numerical results.In this paper, two new interpolation functions are proposed and examined; they are proven to be generally applicable and satisfactory.
Key words: Dual-reciprocity boundary-element method, interpolation function.
1 INTRODUCTION
terms of the power series of a distance function, i.e.: (3O
E
The dual-reciprocity boundary-element method (DRBEM) was first proposed by Nardini and Brebbia 1,2 in 1982 and was shown by Partridge and Brebbia 3 in 1990 to be a general approach to solve elliptic differential equations of the type V2u = b, which represent various kinds of governing equations arising from engineering and applied mathematical problems. 4-7 Hence the success of DRBEM is of great significance to engineers and applied mathematicians. In DRBEM, a series of interpolation functions is used to approximate b in order to convert the associated domain integrals appearing in traditional boundaryelement analyses into boundary integrals only. Hence the choice of interpolation functions has direct effects on the accuracy of this method. One possible choice, as used by Partridge and Brebbia, 3 is of the first several
t)n(x' y)'
(1)
n=0
where l)(x, y) denotes the distance between a source point j and a field point (x, y). According to Partridge and Brebbia, 3 the adoption of a comparatively simple form of interpolation function, i.e. the first two terms of (1) (1 + ~), gives the best results. However, as shown earlier by the present authors, s if b is a function containing partial-derivative terms of the unknown function u(x, y), such a choice results in the creation of singularities on all boundary and internal nodes used in a dual-reciprocity boundary-element analysis. Further evidence is provided in this paper to demonstrate such uncertainty and unreliability associated with the DRBEM. To avoid these artificially created singularities, Zhu and Zhang s proposed a transformation that 387
388
Y. Zhang, S. Zhu
eliminates all derivative terms in b and thus obtained much improved numerical results. Unfortunately, this transformation is not applicable to non-linear governing equations. In this paper, we shall show that, by slightly modifying the originally proposed interpolation functions, not only can much improved results be obtained for the case when b is a function containing partialderivative terms of the unknown function u(x, y), but also the same level of numerical accuracy, if not much better, can also be achieved for other cases in which b does not contain partial-derivative terms of the unknown function u(x,y). The new interpolation functions are therefore highly recommended simply because of their reliability for all the cases. In the next section, the mathematical formulation of the DRBEM is briefly outlined before an example is given to illustrate the large numerical errors resulting from the artificially created singularities introduced by the interpolation functions used by Partridge and Brebbia. 3 One way to avoid the artificially created singularities is suggested with the second term in eqn (1) being dropped so that the interpolation functions are of second-order differentiability. Then, in Section 3, the numerical results obtained upon adopting the two new forms of interpolation function are compared with those obtained from adopting the traditional interpolation functions for all the cases discussed by Partridge and Brebbia. 3 The purpose of showing such a comparison, as clearly stated in the previous paragraph is, of course, not to try to demonstrate that the new interpolation functions can generate much more accurate results than the traditional interpolation functions for all the cases discussed by Partridge and Brebbia 3 as well. On the contrary, our conclusion, which will be given in the last section, is that the new interpolation functions may be preferable for all the cases simply because, for the cases in which singularities are likely to be generated, the new interpolation functions can lead to much improved numerical results, whereas, for the cases in which there is no problem in adopting the traditional interpolation functions, at least the same level of numerical accuracy, if not better, can be achieved by adopting the new interpolation functions.
2 FUNDAMENTALS OF DRBEM AND TWO NEW C H O I C E S OF I N T E R P O L A T I O N FUNCTION The fundamental mathematical theories of DRBEM have been explained in detail by Partridge and Brebbia. 3 Only its main points are therefore outlined here. Consider the following general case of an elliptic governing equation: V2u = b in f~,
(2)
where V 2 is the Laplacian operator, b may be a function
Fig. 1. An illustration of a computational domain and collocation points. of x,y,t,U, Ux, Uy,Ut, and ft is the domain to be considered (see Fig. 1). The DRBEM representation of eqn (2) is given by: 3 Hu -
= (I-rd -
(3)
where H and G are of conventional sense, and each column of Q or (~ contains a vector 6: or cij, respectively. The vector a is calculated via:
a = F-]b,
(4)
where each column of the matrix F contains a vector f:. with its components equal to the values of f:, the interpolation function used in D R B E M , on the corresponding nodes. So far, interpolation functions of the following form are commonly adopted: P
)~ = E
§."(x, y).
(5)
n=O
Finally, substituting eqn (4) into eqn (3) produces the linear system of equations for the unknown u's and q's: Hu - Gq = (HI)
- G(
)F-lb.
(6)
As was shown in detail earlier, 3 the vector a in eqn (3) is either known (in the case when b is a known function of x, y only) or can be associated with the vector u. Thus the unknown u's and q's can be solved from the final linear system of equations. However, if b contains derivatives of u(x, y), the association of Ou/Ox or Ou/Oy to the vector u requires the introduction of a matrix OF/Ox or OF/Oy, since Ou OF,,-1 ~ x x ~ x x r u,
Ou
OF
-1
u.
(7) (8)
As was pointed out by Zhu and Zhang, 2 the diagonal entries of the matrix OF/Ox or OF/Oy are actually undefined. This problem originates from the nondifferentiabifity of the distance function rj(x, y), i.e. the second term in eqn (5), at its source point j. Take the formulation of the matrix OF/Ox as an example. The ith diagonal element of the matrix OF/Ox is defined as
Choice of interpolationfunctions in DRBEM eii= [O~/OX]e~al~t~at point ;, where S = 1 + ri (as suggested by Partridge and Brebbial). But, since Ofi/Ox = 0(1 +ri)/Ox = ( x - xt)/ri and the limit of (x-xi)/ri does not exist as ( x , y ) ~ (xi, yi), eii is actually undefined. Hence expressing Ou/Ox as (OF/Ox)F-1u, which is necessar~ when b contains derivatives of u(x, y) (see other work°'S), creates singularities on all collocation points. To show the numerical inaccuracy resulted from these artificially created singularities, a simple example is given here. Consider the following simple boundaryvalue problem:
50u
v e u = 2 0 x inside a square slab,
(9)
u = e 2~ cosh y on the boundary (see Fig. 2), (10) the exact solution being u = e 2x cosh y. In this case, the vector a in eqns (3) and (4) can be written as: 5 F _ 10u
(11)
where each component of the vector Ou/Ox is equal to the nodal value of Ou/Ox. After substituting eqn (11) into eqn (3) and taking eqn (7) into consideration, the system of equations to be solved for unknown u's and q's can be obtained as: Hu - Gq = ~ (HI] - G(~)F-1 0F F-lu.
(12)
ox
For simplicity, suppose that m - 1 by m - 1 internal nodes (which are required by DRBEM) are evenly distributed in the domain as shown in Fig. 2. Now if Partridge and Brebbia's interpolation functions 3 (referred to as the traditional approach from now on) are adopted, which take the form = 1 + i),
(13)
where J) stands for the distance between the source point
Y 1.0 (
0
1.0 •
Collocationnodes
0
Extreme points of boundaq/element
x
Fig. 2. A sketch of a computational domain (m = 3).
389
j and a field point (x, y), all the diagonal elements of the matrix OF/Ox in eqn (12) are undefined owing to the singular behaviour of the derivative of the interpolation functions near source points. However, if such a nondifferentiability at source points were ignored and all of the diagonal entries were placed zero instead of as shown by Partridge and Brebbia, 3 large numerical errors would be generated when the vector (OF/Ox)F-lu was used to approximate Ou/Ox. This can be clearly seen from the second and fourth columns of Table 1, in which the values of Ou/Ox and its approximating values (OF/Ox)F-lu at 36 collocation points (m = 5) are compared. The approximation of Ou/Ox by (OF/Ox)F-lu leads to a maximum relative error of 70%. In fact, large errors do not only occur at a small number of collocation points; the relative errors on fourteen internal collocation points exceed 20%, and the relative errors are even larger on the boundary collocation points. On the other hand, by using one of the new interpolation functions, which will be presented later (see eqn (15)), the maximum error is just around 7"5%, as shown in the last column in Table 1. Furthermore and most importantly, the numerical errors associated with adopting the traditional approach did not seem to decrease as the number of collocation points was increased, whereas such a reduction was clearly present when the new interpolation function was adopted. As an example, when we increased the number of collocation points from 36 to 121 (m = 10), the maximum error associated with the traditional approach was still about 68 %, and the errors on 23 collocation points exceeded 20%. On the other hand, if the new interpolation function was adopted, the maximum error was further reduced to just 3-5%. It is very clear that the reason for this large discrepancy is the non-differentiability of the interpolation functionfj = 1 + ~) at the source point j, resulting in the diagonal entries of the matrix OF/Ox, undefined since an improved numerical accuracy from 70% to 7"5% was achieved by deleting only one term in the traditional interpolation functions, i.e. the term that leads to singularities in OF/Ox. Furthermore, whether or not the errors from the introduction of artificially created singularities can cancel each other out in the subsequent numerical calculations in solving the boundary-value problem of eqns (9) and (10) was examined. We first solved the above problem by placing linear elements on the boundary. Considerable errors in the calculated flux on the boundary were found with the traditional approach. It was then speculated that the errors might have resulted from the difficulty of representing normals at corner points with finear elements. Linear elements were thus replaced by constant elements. The maximum relative error found for the numerically calculated flux on the boundary, with the traditional approach adopted, was still 46% for m = 5, and it hardly
Y. Zhang, S. Zhu
390
Table 1. Approximation for the vector Ou/Ox
ith node
f2 = 1 + ~*
~ = 1 + r22 + t)3.*
Exact***
Errorlt(%)
Error2:~(%)
1 2 3 4 5 6 7 8 9 10 ll 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
0"180E + 01 0"371E + 01 0"549E + 01 0"833E + 01 0"106E + 02 0"997E + 01 0"108E + 02 0"115E + 02 0"117E + 02 0"102E + 02 0"140E 4- 02 0"128E 4- 02 0'841E 4- 01 0"569E 4- 01 0"289E 4- 01 0"989E + 00 0"158E 4- 01 0"155E + 01 0"132E 4- 01 0"595E + 00 0"321E + 01 0.459E 4- 01 0.690E 4- 01 0"109E 4- 02 0"336E + 01 0.485E + 01 0.726E 4- 01 0.116E + 02 0.368E 4- 01 0-532E + 01 0-798E 4- 01 0.128E + 02 0.422E 4- 01 0.606E + 01 0-915E 4- 01 0"148E + 02
0"245E + 01 0"364E + 01 0"545E + 01 0"807E + 01 0"123E + 02 0"148E + 02 0"152E + 02 0"163E + 02 0"179E + 02 0"198E 4- 02 0"183E 4- 02 0"126E 4- 02 0"837E 4- 01 0"567E 4- 01 0"361E 4- 01 0"280E 4- 01 0"267E + 01 0-241E 4- 01 0"225E 4- 01 0"215E 4- 01 0.300E + 01 0.456E + 01 0.676E + 01 0.101E 4- 02 0.318E + 01 0.484E + 01 0.715E 4- 01 0.108E 4- 02 0.349E 4- 01 0.531E + 01 0.781E 4- 01 0.119E 4- 02 0.397E 4- 01 0.598E + 01 0.878E 4- 01 0.136E 4- 02
0"244E + 01 0"364E + 01 0"544E + 01 0"811E + 01 0"121E + 02 0"149E + 02 0"154E + 02 0"167E + 02 0"185E + 02 0"212E + 02 0'187E 4- 02 0' 125E + 02 0"839E 4- 01 0"562E 4- 01 0"377E 4- 01 0"287E 4- 01 0"251E 4- 01 0"226E + 01 0.209E 4- 01 0.201E + 01 0.304E + 01 0.454E 4- 01 0.677E + 01 0.101E 4- 02 0"323E + 01 0.481E + 01 0-718E 4- 01 0.107E + 02 0.354E 4- 01 0.528E + 01 0.787E 4- 01 0.117E 4- 02 0.399E + 01 0.595E + 01 0-888E 4- 01 0.132E + 02
26 1"9 l'0 2"8 13 33 30 31 37 52 25 24 0"24 1"2 23 66 37 31 37 70 5'4 1.1 1.8 8.4 4.0 0.72 1.2 8"1 4.1 0.77 1"3 9.0 5.7 1.8 3"0 12
0"21 0"017 0"17 0"50 1"3 0"091 1"6 2"5 3"8 6"7 1"9 0"44 0"19 0"75 4"2 2"2 6"2 6"8 7"4 7"0 1"6 0"32 0"17 0"16 1"5 0"51 0"43 0"77 1"4 0"61 0"74 l'l 0"58 0 "46 1'1 2"5
*Numerically calculated Ou/Ox by adopting the interpolation functions proposed by Partridge et al. **Numerically calculated Ou/Ox by adopting the interpolation functions proposed by us. ***Exact value of Ou/Ox. tErrorl = I(Column2- Column4)/Column41 • 100. ~:Error2 = I(Column3- Column4)/Column41 • 100.
declined with m being increased. This led us to conjecture that the artificially created singularities were responsible for the large errors present in the final D R B E M solution. To prove this conjecture, we introduce two new interpolation functions. Deleting the second term in eqn (1), which yields an interpolation function differentiable up to the second order at any point, is obviously a natural solution to avoid creating these singularities. The simplest choice can be either .~ = 1 + t)?
(14)
fj=
(15)
or 1 +~2 +~3.
The expression ~ = 1 + ~2 cannot be adopted because we found that the matrix F becomes non-invertible in
many situations with this form of interpolation function being adopted. U p o n simply replacing eqn (13) by either eqn (14) or eqn (15) while keeping everything else exactly the same, much more accurate numerical results were obtained for both the unknown, u, and the flux, q, even with a smaller number o f internal nodes and that o f boundary elements being used in the previous example. For simplicity, only the results from adopting eqns (13) and (15) as the interpolation functions are presented here. For m = 5, the maximum errors for the numerical solutions on all internal nodes are small (0.9% for using eqn (13) as the interpolation functions but only 0.06% for using eqn (15) as the interpolation functions). On the boundaries, as shown in Table 2, the maximum relative errors for the numerically calculated flux are 46% and 6% for using eqns (13) and (15) as the interpolation functions, respectively (we excluded those boundary nodes for
Choice of interpolation functions in D R B E M
391
Table 2. Numerically calculated and exact flux
Boundary nodes 1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
P-Approach*
N-Approach**
Exact***
ErrorH'(%)
Error2~(%)
-0"305 -0"002323 0"119 0"190 -0"252 14"01 14"95 16"06 17"58 20"20 6"645 4"873 3"364 2"203 0"7688 -2"570 -2"712 -2"358 -2"306 - 1"572
0"0122 0"00866 0"00275 -0"00869 0"0345 14"85 15"43 16"64 18"48 21 "28 7"051 4"770 3"172 2"130 1"451 -2"806 -2"504 -2"239 -2"056 -2" 130
0" 0" 0" 0" 0" 14"85 15"45 16"66 18"55 21"18 7"110 4"766 3"195 2"141 1'435 -2"866 -2"510 -2"255 -2"091 -2"010
5"7 3"3 3"6 5"2 4"6 6"5 2"3 5"23 2"9 46 10 8"0 4"6 10 21
0"026 0"12 0"17 0"39 0"50 0"82 0"094 0"72 0"52 1"1 2"1 0"25 0"72 1"7 6"0
*Numerically calculated flux by using the traditional approach proposed by Partridge et al. **Numerically calculated flux by using our new approach. ***Exact solution. tErrorl = I(Column2 - Column4)/Column41 • 100. ~Error2 = I(Column3 - Column4)/Column41 * 100.
which the exact solutions are equal to 0). Once again, when the total number of coUocation nodes and that of the boundary elements were increased, the numerical errors associated with adopting eqn (13) as the interpolation functions hardly declined. For example, when m = 10, there were still three points at which the numerical errors exceeded 20*/0 (one of them even up to 57%!), not even counting the errors on those boundary points at which the exact solutions are equal to 0. Thus, the only explanation for such a considerable reduction obtained by just dropping one term in the interpolation functions and keeping all the others the same is that the artificially created singularities with eqn (13) being used as the interpolation functions are responsible for the large numerical errors present both in approximating Ou/Ox by (OF/Ox)F-lu and in the final DRBEM solution.
3 N U M E R I C A L EXPERIMENTS AND DISCUSSION With the example shown in the previous section and the other two examples shown elsewhere by the present authors, s we believe that the simplest and most preferable replacements to the interpolation function used by Partridge and Brebbia, 3 i.e. 1 + ~, will be those shown in eqns (14) and (15) when b contains derivatives of u(x, y). For other cases, eqn (13) remains the simplest one and can still be used to generate quite accurate numerical results. However, in order to study the applicability of eqns (14) and (15) to other cases as
well, in addition to adding a new case in Subsection 3.1, we followed Partridge and Brebbia 3 to examine the numerical accuracy associated with adopting eqns (14) and (15) as interpolation functions in those cases presented by them. The results obtained from applying the interpolation functions used by Partridge and Brebbia, 3 the functional transformation suggested by Zhu and Zhang, 8 and the two new interpolation functions proposed in this paper are compared in this section. In the following discussion and in all the tables presented, for easy reference, we define the following approaches: (i) approach 1: the approach with J} being of the form suggested by Partridge and Brebbia; 3 (ii) approach 2: the approach with a functional transformation suggested by Zhu and Zhang; 8 (iii) approach 3: the approach with ~ being of the form shown in eqn (14); (iv) approach 4: the approach with f2 being of the form shown in eqn (15). To compare the four approaches with each other more clearly, we need some measurements to judge the overall performance of each approach. It seems reasonable to adopt as such measurements the maximum absolute error between two functions f and g (they may represent the exact and numerical solutions, respectively), which is defined as: IIf, g II = max
(x,y) • fl
If(x, y) - g(x, y) I,
(16)
Y. Zhang, S. Zhu
392
Table 3. Numerical selatlem for the Pobmn-type problem with approaches 3 and 4
ui*
x
y
Approach 3
Approach 4
Exact solution
Absolute difference 1'*
Absolute difference 2***
u17 ul8 Ul9
1"5 1.2 0.6 0"0 -0.6 -1.2 - - 1"5 -1"2 -0"6 0"0 0"6 1'2 0"9 0"3 0"0 -0'3 -0"9
0"0 -0.35 :0.45 -0.45 -0.45 -0'35 0"0 0"35 0"45 0"45 0"45 0"35 0"0 0"0 0"0 0"0 0"0
1'2954 1.0916 1.0018 1.0242 1.1393 1.3131 1"2954 1"0916 1"0018 1'0242 1'1393 1"3131 1"1018 l'OlO0 0"9987 l'OlO0 l" 1018
1"2952 1.0913 1.0019 1.0244 1.1397 1-3139 1"2952 1"0912 1"0019 1"0244 1"1397 1"3139 1"1020 1"0103 0"9990 1"0103 l" 1020
1"2947 1.0917 1.0028 1.0254 1.1410 1-3156 1'2947 1'0917 1'0028 1'0254 1"1410 1"3156 l" 1030 l'O113 l'O000 1"0113 l" 1030
0.0007 0.0001 0.0010 0.0012 0.0017 0.0025 0.0007 0.0001 0-0010 0"0012 0"0017 0"0025 0"0012 0"0013 0"0013 0"0013 0"0012
0.0005 0.0004 0.0009 0.0010 0.0013 0.0017 0"0005 0"0005 0"0009 0"0010 0"0013 0'0017 0'0010 0"0010 0"0010 0"0010 0"0010
Uzo u21 u22 u23 u24 u25 U26 U27 u28 u29 u30 ual u32 u33
*The values of u(x, y) on the ith node. **Differences between the function values computed by using approach 3 and those from the exact solution. ***Differences between the function values computed by using approach 4 and those from the exact solution.
and the m a x i m u m relative error, which is defined as: max
] (g(x,
y) - f ( x , y))/f(x, y)[,
(f(x, y) ¢ 0),
( x , y ) E f~
(17) where f~ is the domain within which both functions are defined. It is evident that the smaller the m a x i m u m relative error and the m a x i m u m absolute error are, the better an approach is. In the following numerical experiments, the domains are all taken to be an ellipse with semi-major axis of length 2 and semi-minor axis of length 1, as shown in Fig. 3. The boundary is discretized into 16 linear elements. For ease of comparison, the numerical solutions resulting from adopting the two new approaches (i.e. approaches 3 and 4) and the m a x i m u m absolute errors and m a x i m u m relative errors resulting from the four approaches are all tabulated in Tables 3-13.
Y 1.0
•
24Q !3
27
25 33
28
•
3.1 Poisson-type problem In this case, b is a known function of x, y only. The linear system for u's and q's is just eqn (6) (note that the right-hand side of eqn (6) is a known vector). Since b is not a function of Ux and Uy, approach 2 is not applicable. The numerical solutions at the internal nodes with the two new approaches for a particular case b = [cosh{(x+y)/2}]/2 (with the exact solution to the differential equation (2) being u(x, y) = cosh (x + y)/2) are tabulated in Table 3. F r o m the last two columns of Table 3, we can clearly see the good performance of these two new approaches. In this case, approach 4 seems better than approach 3. To examine the overall performance of each approach, the m a x i m u m absolute errors and m a x i m u m relative errors o f approaches 1, 3, and 4 are shown in Table 4, from which we can clearly see the superiority of the two new approaches over the traditional one; the m a x i m u m absolute errors of the two new approaches are twice or four times smaller than that of approach 1 respectively. Similarly, the improved m a x i m u m relative errors of the two new approaches are also twice or four times smaller than that of approach 1, respectively. All these indicate that the functional distance between the numerical solutions from the two new approaches and the exact solution is much smaller.
32 ;'0
O19
Fig. 3. An elliptical domain and a discretization.
Table 4. Overall performance of each approach for the Poissontype problem Norm Maximum relative error
Approach 1
Approach 3
Approach 4
0-0205
0-00899
0.00571
1.6%
0.71%
0.45%
Choice of interpolation functions in DRBEM
393
Table 5. Numerical solutions for the case b = - u with approaches 3 and 4
ui*
x
y
Approach 3
Approach 4
Exact solution
Absolute difference 1"*
Absolute difference 2"**
u17 uls
1.5 1.2
0.0 -0.35
0.9916 0-9318
0-9904 0.9318
0.9975 0.9320
0.0059 0.0002
0.0071 0.0002
Z/19
0"6
--0"45
0"5647
0"5648
0"5646
0"0001
0"0002
u20 u29 u30 U31
0"0 0"9 0"3 0"0
-0"45 0"0 0"0 0"0
0-0000 0-7827 0"2955 0"0000
0.0000 0"7826 0"2955 0"0000
0.0000 0"7833 0"2955 1"0000
0"0000 0"0006 0.0000 0"0000
0"0000 0"0007 0"0000 0"0000
*The values of u(x, y) on the ith node. **Differences between the function values computed by using approach 3 and those from the exact solution. ***Differences between the function values computed by using approach 4 and those from the exact solution.
3.2 b =
b(u), linear case
In the second numerical experiment, b is a function of the unknown function u, but not a function of its derivatives. A particular exact solution for V2u = - u is u(x,y) = sin x, which was used by Partridge and Brebbia 3 to verify D R B E M . Now, the fight-hand side of eqn (6) is simply - ( I l l ] - G{~)F-lu, which introduces no new unknowns but the unknown u's that have appeared on the left-hand side o f eqn (6), so unknown u's and q's can be obtained from eqn (6). The numerical solutions at internal nodes with the two new approaches are shown in Table 5 (owing to the symmetry of the problem, only the solutions on seven internal nodes are shown in Table 5). A good agreement between the numerical solutions and exact solutions is clearly indicated by the last two columns of Table 5. In this case, approach 3 seems to be better than approach 4. This is confirmed by Table 6; the maximum absolute error and maximum relative error o f approach 3 are both smaller than those of approach 4. In this case, the maximum absolute errors of approaches 1, 3 and 4 (see Table 6) are all very small and close to each other, as are the maximum relative errors. Again, good numerical accuracy resulting from adopting the two new approaches has been demonstrated for this case. 3.3 b =
b(x, y, u, ux, uy),
finear c a s e
As was stated in Section 2, if b contains derivatives of u(x, y), approach 1 becomes unreliable. If b is still linear in terms of the unknown function u and its derivatives, approach 2 can be employed to avoid artificially created singularities.
Table 6. Overall performance o f each approach for the case
b = -u Norm Maximum relative error
Approach 1
Approach 3
Approach 4
0"0107
0.0131
0.0159
1.7%
1.3%
1.6%
Let us examine two linear cases first: b = - O u / O x , with the exact solution u(x, y) = e -x, and b = x(Ou/Ox), with the exact solution u(x, y) = xe y (for the latter case, if the example used in Subsection 4.3 of Ref. 3, i.e. b = - ( 2 / x ) (Ou/Ox), is adopted, the transformed governing equation with approach 2 is a Laplace's equation, which is a trivial case; the example used in Ref. 8 is therefore adopted here). With approaches 1, 3, and 4, it can be shown 3,8 that the right-hand side of eqn (6) can be written in terms of the vector u as - ( H I ] - G O ) F _ 1 ~xxrOF- l u and
(HIll -- GO)F-1X ~ - F - l u , respectively. In these expressions, the matrix
Ox =
L
(xi, y,)
(with approach 1, the diagonal entries o f the matrix OF/Ox, which are actually undefined, s are artificially set to be zero, as was done in Ref. 3), and X is a diagonal matrix with its non-zero entries being equal to the value of the x-coordinate at each node. With approach 2 , a transformation u = r e - x/2 or 1./= ~'ex2/4, where v is a new unknown function, can be constructed to a v o i d dealing with the derivative terms directly. The vector b on the right-hand side o f eqn (6) can then be written directly in terms of the vector v. The numerical solutions at internal nodes with the two new approaches are shown in Tables 7 and 8, respectively. Again, the last two columns of Tables 7 and 8 show the accuracy of the numerical solutions resulting from adopting the two new approaches. The maximum absolute errors and maximum relative errors of each approach are given in Tables 9 and 10. For the case b = -Ou/Ox, it seems that approach 2 produces the best numerical results and that the two new approaches are better than the traditional approach. For the case b = x(Ou/Ox), the two new approaches give much smaller maximum absolute errors than the other two approaches, which again shows the
Y. Zhang, S. Zhu
394
Table 7. N ~ m r l e a l solutions for the ease b = - ~ / S x
with ~
3 and 4
ui*
x
y
Approach 3
Approach 4
Exact solution
Absolute difference 1"*
Absolute difference 2***
U!7 Ula ul9 u20 u21 u22 U23 1/29 U30 U31 U32 U33
1"5 1"2 0.6 0-0 -0"6 - 1.2 1-5 0"9 0"3 0"0 --0"3 --0"9
0"0 --0"35 -0.45 -0"45 -0"45 -0.35 0"0 0"0 0"0 0"0 0"0 0"0
0"2253 0"2994 0.5480 1.0017 1.8284 3.3297 4.4779 0"4056 0"7411 1"0015 1"3528 2"4663
0"2265 0"2987 0"5469 1.0006 1"8273 3"3286 4-4779 0"4050 0"7401 1"0005 1"3517 2"4654
0"2231 0"3012 0"5488 1.0000 1"8221 3.3201 4.4817 0"4066 0"7408 1"0000 1"3498 2"4596
0"0022 0"0018 0.0008 0.0017 0.0063 0.0096 0.0038 0"0010 0"0003 0"0015 0"0030 0"0067
0"0034 0"0025 0.0019 0.0006 0.0052 0.0085 0.0038 0"0016 0"0007 0"0005 0"0019 0"0058
-
-
*The values of u(x, y) on the ith node. **Differences between the function values computed by using approach 3 and those from the exact solution. ***Differences between the function values computed by using approach 4 and those from the exact solution.
Table 8. Numerieal solutions for the ease b = x ( ~ / S x ) with approaches 3 and 4 ui*
x
y
Approach 3
Approach 4
Exact solution
Absolute difference 1"*
Absolute difference 2***
ul7 uls ux9 u20 U26 U27 u28 u29 U3o u31
1-5 1"2 0"6 0"0 --0"6 1"2 - 1"5 0"9 0"3 0-0
0-0 -0"35 -0"45 -0"45 --0"45 -0"35 0"0 0"0 0"0 0-0
1"4987 0"8355 0"3764 0.0000 0"0000 0"9268 1"6802 0"8888 0-2959 0.0000
1"4992 0"8326 0"3757 0.0000 0"0000 0-9274 1-6827 0"8888 0"2959 0.0000
1"5000 0"8456 0"3828 0.0000 0"0000 0"9410 1"7029 0"9000 0"3000 0.0000
0"0013 0"0101 0-0064 0.0000 0"0000 0"0142 0"0227 0"0112 0.0041 0.0000
0"0008 0"0130 0.0071 0.0000 0"0000 0"0136 0"0202 0.0112 0"0041 0.0000
-
-
*The values of u(x, y) on the ith node. **Differences between the function values computed by using approach 3 and those from the exact solution. ***Differences between the function values computed by using approach 4 and those from the exact solution.
Table 9. Overall performnee of each mEgomch for the case b =
Norm Maximum relative error
Approach 1
Approach 2
Approach 3
Approach 4
0-0383 7.1%
0.0134 0.4%
0.0269 3.6%
0.0247 5-4%
Table 10. Overall performmmce of each mEpnm~h for the c ~ e b =
Norm Maximum relative error
-Ou/Ox
x(Su/Sx)
Approach 1
Approach 2
Approach 3
Approach 4
1-03 6.7%
1.37 3.7%
0-249 4-3%
0.383 5.3%
Choice of interpolation functions in DRBEM
395
Table 11. Nmme~al solultons for the ease b = - u ( ~ / S x ) with apFm~m~es3 and 4 u~*
x
y
Approach 3
Approach 4
Exact solution
Absolute difference 1'*
Absolute difference 2***
ul7 uls u19 u20 u21 u~ u23 u2o u30 Ual u32 u33
1"5 1.2 0.6 0-0 -0"6 - 1.2 - 1"5 0"9 0"3 0.0 -0.3 -0.9
0"0 -0.35 -0.45 -0.45 -0.45 -0.35 0.0 0.0 0.0 0-0 0.0 0"0
0.4531 0.4716 0.5498 0.6629 0.8329 1.1145 1.3309 0.5089 0.6016 0"6630 0"7382 0.9533
0.4566 0.4707 0-5476 0.6606 0.8307 1.1127 1.3307 0.5076 0"5995 0.6609 0.7361 0-9515
0"A.A.A.A. 0.4762 0.5555 0.6667 0.8333 1.1111 1-3333 0-5128 0-6061 0"6667 0.7407 0.9524
0.0087 0.0046 0.0057 0.0038 0.0004 0.0034 0.0024 0.0039 0.0045 0.0037 0.0025 0-0009
0.0122 0.0061 0.0079 0.0061 0.0026 0.0016 0.0026 0.0052 0.0066 0.0058 0.0046 0.0009
*The values of u(x, y) on the ith node. **Differences between the function values computed by using approach 3 and those from the exact solution. "**Differences between the function values computed by using approach 4 and those from the exact solution. high degree o f accuracy o f the two new approaches. Furthermore, in all the numerical experiments we have conducted, it was found that the two new approaches produced consistent and accurate results, whereas approach 1 may sometimes yield very large numerical errors, as shown in the example given in the previous section. Hence, for the case in which b contains derivatives o f u(x, y), the two new forms o f interpolation function are much more reliable choices for DRBEM.
3.4 b = b ( x , y , u . , uy), non-rmear case
The final numerical experiment is a non-linear case, with b = -u(Ou/Ox) and the exact solution being u(x, y) = 2/x. F o r this case, approach 2 does not apply. To avoid the singularity in the solution on some nodes, a change o f the origin in Fig. 3 to ( - 3 , 0) was made. 3 In this case, the right-hand side of eqn (6) can be written as
- ( t I C - C,Q)r-
v
where the diagonal o f the matrix U contains the nodal values o f u(x, y). Since this is a non-linear case, an iterative procedure has to be used to solve for unknown u's and q's. Initially, the Laplace equation H a - Gq = 0 is solved to produce a Ul vector, which is put back to the diagonal of U; then the equation system eqn (6) is solved to generate u2, and so on, until a convergent result is obtained. The numerical
Table 12. Overall performnce of each approach for the ease b = -,,(o,,/ax)
Norm Maximum relative error
Approach 1
Approach 3
Approach 4
0.0190
0.0216
0-0304
2.4%
5.1%
7-1%
solutions on twelve internal nodes, with approaches 3 and 4 being adopted and with the convergence criterion being set to be 1% (i.e. when the change for the calculated u's and q's is less than 1%, non-linear iterations are stopped), are shown in Table 11. The results are still quite satisfactory. Table 12 shows the maximum absolute errors and maximum relative errors for each approach. Approach 1 seems to have produced the best results. However, further examination revealed that the calculated q's on the boundary were of very large errors (some o f them even more than 100%D and that the numerical results from this particular nodal pattern (see Fig. 3) were unstable, i.e. when some minor changes were given to the nodal pattern, a considerable amount o f change took place in the numerical solutions. More boundary elements and internal nodes were therefore included until stable results had been obtained. The new maximum absolute errors and maximum relative errors for each approach are tabulated in Table 13, from which we can clearly see that the results obtained from adopting approaches 3 and 4 are a little more accurate than those from adopting approach 1.
4. C O N C L U S I O N In this paper, it is dearly shown that the artificially created singularities with the traditional interpolation functions being used in the dual-reciprocity boundaryelement method can result in large numerical errors in
Table 13. New norms and maximum relative errors for the case b = -,,(a./ax)
Norm Maximum relative error
Approach 1
Approach 3
Approach 4
0-00867
0.00277
0.00334
0.65%
0.46%
0.59%
396
Y. Zhang, S. Zhu
the case in which b contains derivatives of the unknown. Although large numerical errors may sometimes not be present in the final numerical solution (e.g. in the particular non-linear example shown in subsection 3.4, reasonably accurate numerical results were still obtained), D R B E M with the adoption of the traditional interpolation functions becomes unreliable, since several examples with large numerical errors have been reported by Zhu and Zhang 8 and are shown in the present paper. Two new forms of generally applicable interpolation functions have been proposed to be the best and simplest replacements for the traditional ones when large errors may be generated from artificially created singularities. The results of our numerical experiments suggest that the numerical errors associated with the artificially created singularities can be avoided; much more reliable numerical results have been presented for both linear and non-linear cases. Furthermore, the newly proposed interpolation functions have also been tested by applying them to all the other cases used by Partridge and Brebbia 3 when they presented DRBEM as a powerful solver for elliptic equations of the type V2u = b. Very satisfactory results were also obtained, which does seem to be in favour of replacing the traditional interpolation functions by the newly proposed ones when D R B E M is adopted.
2.
3.
4.
5.
6.
7.
8.
REFERENCES 9. 1. Nardini, D. & Brebbia, C. A. A new approach to free vibration analysis using boundary elements. In Boundary
Element Methods in Engineering, ed. C. A. Brebbia, Springer-Verlag, Berlin, Germany, and New York, NY, USA, 1982. Nardini, D. & Brebbia, C. A. Boundary integral formulation of mass matrices for dynamic analysis. In Topics in Boundary Element Research, Vol. 2, Springer-Verlag, Berlin, Germany, and New York, NY, USA, 1985, Chapter 7. Partridge, P. W. & Brebbia, C. A. Computational implementation of the BEM dual reciprocity method for the solution of general field equations. Communications in Applied Numerical Methods, 1990, 6, 83-92. Loeitter, C. F. & Mansur, W. J. Dual reciprocity boundary element formulation for potential problems in infinite domains. In Boundary Elements X, Vol. 2, Computational Mechanics Publi~tions, Southampton, 1988. Wrobel, L. C., Telles, J. C. F. & Brebbia, C. A. A dual reciprocity boundary element formulation for axisymmetric diffusion problems. In Proceedings of VIII International Conference on BEM in Engineering, Tokyo, Computational Mechanics Publications, Southampton, UK, and Boston, MA, USA, 1986. Brebbia, C. A. & Wrobel, L. C. The dual reciprocity boundary element formulation for non-linear diffusion problems. Computational Methods in Applied Mechanics and Engineering, 1987, 65, 147-64. Partridge, P. W. & Brebbia, C. A. The dual reciprocity boundary element method for the Helmholtz equation. In Mechanical and Electrical Engineering, ed. C. A. Brebbia & A. Chaudouet-Miranda. Computational Mechanics Publications, Southampton, UK, and Boston, MA, USA, 1990. Zhu, S. & Zhang, Y. An improvement on the dual reciprocity boundary element method for equations with convective terms. Communications in Applied Numerical Methods, in press. Brebbia, C. A., Telles, J. C. F. & Wrobel, L. C. Boundary Element Techniques, Springer-Verlag, Berlin, Germany and New York, NY, USA, 1984.