- Email: [email protected]
A modified Gauss quadrature formula with special integration points for evaluation of Quasi-singular integrals
ARTICLE IN PRESS
Engineering Analysis with Boundary Elements 30 (2006) 758–766 www.elsevier.com/locate/enganabound
A modified Gauss quadrature formula with special integration points for evaluation of Quasi-singular integrals E.E. Theotokoglou, G. Tsamasphyros Faculty of Applied Sciences, Department of Mechanics, Laboratory of Strength of Materials, The National Technical Universtiy of Athens, Zoghraphou Campus, Theocaris Bld., GR-0157 73, Athens, Greece Received 11 November 2005; accepted 13 May 2006 Available online 14 July 2006
Abstract It is well known in the boundary element method that integration rules fail when the integrand presents a nearby singularity. This drawback arises when the field point is near the source point, i.e. in the case of a domain with very narrow boundaries or when the field point where we try to calculate stresses or any other field variables, is near the boundaries. In the present paper a quadrature formulas for ‘ isolated singularities near the integration interval, based on ordinary or special Langrange interpolatory polynomials, is obtained. This interpolatory formulas present similarities with known formulas for the numerical evaluation of singular integrals. Quadrature formulas for regular and singular integrals with conjugate poles are also derived. Numerical examples are given and the proposed quadrature rules present the expected polynomial accuracy. r 2006 Elsevier Ltd. All rights reserved. Keywords: Nearby poles; Quadrature formula; Lagrangian interpolation; Cauchy type integrals; Elasticity; Boundary element method; Numerical integration
1. Introduction In the solution of two- and three-dimensional problems with the boundary element method (BEM) [1–15], the BEM is required to evaluate integrals where the integrand is of the form 0(‘n(r)), 0(1/r) or 0(1/r2) in which r is the distance between the source point and a field point. If the source point belongs to the integration interval we say that the integral has a (i) weak singularity if the kernel behaves like ‘nr or 1=rg with go1 and the integrals are defined in the usual way. (ii) A strong singularity if the kernels behave like 1=rg with gX1. Integrals of this type are diverging integrals and can be defined as Hadamard finite part integrals [4–8]. If g ¼ 2k þ 1, k ¼ 0; 1; 2; . . ., the corresponding integrals can be defined as Cauchy principal value integrals [6–8]. Corresponding author. Sofouli 33, Nea Smirni 17122, Athens, Greece. Tel.: +30 210 772 1303; fax: +30 210 7721302. E-mail address: [email protected] (E.E. Theotokoglou).
0955-7997/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2006.05.001
Up to now many methods are adopted for the calculation of such integrals. So the integral is calculated by splitting the boundary in smaller or bigger integration intervals. The numerical evaluation of such integrals is very crucial in the BEM, not only for singular integrals but also for the regular ones, in order to compute field variables and their derivatives. The kernels of integral forms present a similar behavior ð1=rg Þ giving primary or secondary variables (like displacements, strains and stresses in the case of elasticity). Taking into consideration that secondary variables are mth order derivatives of primary variable entered on the boundary integral equation, their behaviour is like 1=rg with gX1. If the source point does not belong to the integration interval the integral is defined as a regular one. But if the distance r, of the source point and the nearest field point is small, the integrand takes a great value. The difficulties arise from the fact that the integrands of integrals along the boundary vary rapidly with the distance when the field point approaches the source point. If r is very small this error is very important and surprisingly is much more important as the number of integration points increase. We
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
say that these integrals present ‘‘a nearby singularity’’ or a quasi-singularity. In engineering problems, these difficulties arise in the analysis of thin or shell-like bodies [9,10], in bodies with very narrow boundaries, i.e. in I-section beams or in the necking problem or in curved thin members, in cracks near boundaries [11–13], in contact problems [14], and in the case where the boundary integral is constructed in order to calculate displacements and stresses at internal points, near the bounding surface of the body. In the particular case of a formulation based on Muskhelishvili’s complex potentials [1,11,15] in a plane elastic body D bounded by two smooth non-intersecting arcs L1 and L2 ðL ¼ L1 [ L2 Þ, the complex stress vector rn ðz0 Þ at a source point z0, is given by rn ðz0 Þ ¼ snn ðz0 Þ þ isnt ðz0 Þ ¼ Fðz0 Þ þ Fðz0 Þ o d¯z0 n z0 F0 ðz0 Þ þ Cðz0 Þ ; z0 2 D, þ dz0
ð1Þ
where snn ðz0 Þ is the normal component of the stress vector along an arbitrary fictitious cut ðt tÞ, snt ðz0 Þ is the corresponding tangential stress, and Fðz0 Þ, Cðz0 Þ are the Muskhelishvili’s complex potentials. These complex potentials can be expressed in terms of a distribution of complex dislocations g(t) [1,6,11] along L1 and L2, by the following singular integrals Z Z 1 gðtÞ 1 gðtÞ Fðz0 Þ ¼ dt þ dt, (2) 2pi L1 t z0 2pi L2 t z0 Cðz0 Þ ¼
Z Z 1 gðtÞ 1 gðtÞ dt þ dt 2pi L1 t z0 2pi L2 t z0 Z Z ¯tgðtÞ ¯tgðtÞ þ dt þ dt. 2 2 L1 ðt z0 Þ L2 ðt z0 Þ
If z0 is very close to L1 or L2 these integrals behave like 1=r ¼ j1=ðt z0 Þj, i.e. the integrals have a nearby singularity. Let now L1 and L2 are very narrow. If z0 belongs to L1 the integrals in (2) along L2 behave like nearby singular integrals and vice versa. The general form of an integral with nearby and strong singularities can be put in the following form: Z 1 wðtÞjðtÞ dt, (3a) IðjÞ ¼ 1 oðtÞ with oðtÞ ¼
‘ Y ðt yi Þ;
yi ¼ xi þ ii ,
(3b)
759
length s of a part of the boundary or of the whole boundary of the body D, can be transformed to an integral along ð1; 1Þ. Among the numerical methods for evaluation of nearby singularities, regularization procedures [16–19] and procedures trying to avoid singularities [20–24] appear. On the other hand, procedures based on transformations, which remedy the singularities [25–28] and on procedure based on general quadrature rules [29–44] have also been proposed. Strong, weak and nearby singularities have already been confronted at early eighties by our team [34,37,38,41]. Especially nearby singularities are confronted in two master thesis [37,38] and the results are incorporated in our codes for the solution of boundary integral equations. In the present study, we have expanded the known Gauss-integration rules in order to confront the above problem where many strong or weak singularities can coexist. For the sake of generality we have considered the most general case by supposing that, for some reasons, we need to include some preassigned nodes. A ‘‘modified’’ or ‘‘pseudo-Gauss’’ quadrature formulas has been proposed. This quadrature formula is produced from a simple Lagrange interpolation jðtÞ. For question of generality preassigned nodes are considered, so Lobatto, Radau and Kronrod formulas [34,35,44] are derived. Based on the proposed rule, we have also constructed modified quadrature formulas for regular integrals with a single pole, for regular integrals with conjugate poles: Z 1 wðtÞjðtÞ dt; c51, 2 2 1 t þ c where the polynomial accuracy would be equal to 2n+1 for n integration points and for Cauchy singular integral with conjugate poles: Z 1 wðtÞjðtÞ dt; c51; 1oxo1, 2 2 1 ðt xÞðt þ c Þ where the polynomial accuracy would be equal to 2n+2. Thus, a greater polynomial accuracy has been achieved in the case of the proposed formula, in comparison with other methods and procedures [32–35]. Finally, examples are presented for the application of the modified quadrature rule illustrating, regular and singular Cauchy type integrals with conjugate poles. 2. Some definitions
i¼1
where ei is small or zero and xi is on ½1; 1 or very close to this interval and w(t) the weight function which is fixed positive and integrable on ½1; 1, given by wðtÞ ¼ ð1 tÞa ð1 þ tÞb ;
a; b4 1.
(3c)
The previous form (3) concerns not only line integrals but also surface integrals if the latter can be reduced to one-dimensional integrals and then the form (3) can be applied. Obviously any integral with respect to the arc
Let us give some necessary definitions: We denote by D a simply connected domain containing ð1; 1Þ in its interior. We denote by jðzÞ the analytic continuation of jðtÞ into D. Let C a simple contour inside D, presenting angular points at 1 and 1 with angles a and b, respectively (Fig. 1). These angular points are due to the behaviour of jðzÞ.
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
760
or in more compact form. jnþmþ‘ ðtÞ ¼
nþm X j¼1
þ
oðtÞPnþm ðtÞ jðtj Þ ðt tj Þoðtj ÞP0nþm ðtj Þ
‘ X i¼1
oðtÞPnþm ðtÞ jðyi Þ ðt yi Þo0 ðyi ÞPnþm ðyi Þ
ð8Þ
and Fig. 1. A simple connected domain D containing (1, 1), at its interior.
oðtÞPnþm ðtÞ Rnþmþ‘ ðtÞ ¼ 2pi
I
jðzÞ dz. ðz tÞoðzÞPnþm ðzÞ
(9)
C
ftj gnj¼1 ; fzk gm k¼1
Let two sets of points inside D, belonging or not to the integration interval. We suppose that zk are preassigned, whereas tj can be selected arbitrarily. These points will be used as interpolation points. Besides these interpolation points, it is possible to use as ‘‘special integration points’’ the roots yi ði ¼ 1; 2; . . . ; ‘Þ of oðtÞ ¼ 0. In this case the interpolation polynomial jnþmþ‘ ðzÞ is a complex polynomial of degree ðn þ m þ ‘ 1Þ, and we must have jðtj Þ ¼ jnþmþ‘ ðtj Þ; where ( tj ¼
tj ; zk ;
jðyi Þ ¼ jnþmþ‘ ðyi Þ,
(4)
1pjpn; j ¼ n þ kpn þ m:
We denote by pn ðtÞ ¼ Y nþm
n Y
ðt tj Þ;
OðtÞ ¼
j¼1
m Y
ðt zk Þ,
k¼1
ðtÞ ¼ pn ðtÞOðtÞ;
W ðtÞ ¼ wðtÞOðtÞ.
ð5Þ
3. Quadrature formulas with special integration or collocation points (modified or pseudo-Gauss formula with correction terms) As we mentioned above, we may use as special integration points the roots yi ði ¼ 1; 2; . . . ; ‘Þ of oðtÞ. Thus from the interpolation formula (4) and taking into consideration Lagrange procedure [45–47], it is obtained jðtÞ ¼ jnþmþ‘ ðtÞ þ Rnþmþ‘ ðtÞ, jnþmþ‘ ðtÞ ¼
n X j¼1
þ þ
(6)
oðtÞPnþm ðtÞ jðtj Þ ðt tj Þoðtj ÞOðtj Þp0n ðtj Þ
n X
oðtÞPnþm ðtÞ jðzk Þ ðt zk Þoðzk ÞO0 ðzk Þpn ðzk Þ k¼1
‘ X i¼1
oðtÞPnþm ðtÞ jðyi Þ, ðt yi Þo0 ðyi ÞPnþm ðyi Þ
Multiplying both sides of (6) by wðtÞ=oðtÞ and integrating on ½1; 1, we get Z 1 wðtÞ jðtÞdt ¼ QpG (10) nþmþ‘ ðjÞ þ E nþmþ‘ ðjÞ, oðtÞ 1 where QpG nþmþ‘ ðjÞ ¼
lj
j¼1
‘ jðtj Þ X þ ln jðyi Þ oðtj Þ i¼1 i
is the quadrature formula, Z 1 2c ðtj Þ 1 W ðtÞpn ðtÞ lj ¼ 0 dt ¼ 0 n Pnþm ðtj Þ 1 t tj Pnþm ðtj Þ
(11)
(12)
are the classical weights of the standard interpolatory formula for regular integrals, Z 1 1 W ðtÞpn ðtÞ n li ¼ 0 dt o ðyi ÞPnþm ðyi Þ 1 t yi 2cn ðyi Þ ð13Þ ¼ 0 o ðyi ÞPnþm ðyi Þ are the weights corresponding to the special integration points, and Z 1 I 1 jðzÞ W ðtÞpn ðtÞ dt dz E nþmþ‘ ðjÞ ¼ 2pi oðzÞPnþm ðzÞ 1 z t C I 1 jðzÞcn ðzÞ ¼ dz ð14Þ pi oðzÞPnþm ðzÞ C
is the error of the quadrature formula. The first term of (11) can also be obtained if the classical integration formula is applied to the integrand jðtÞ=oðtÞ. If, in other words, the singularities (the poles) of the integrand are ignored, the second term can be considered as a ‘‘correction term’’. The error term following the prescribed procedure and Guelfond [48], can be written as E nþmþ‘ ðjÞ ¼
ð7Þ
nþm X
1 X
mnk1 jðnþmþ‘þkÞ ðrk Þ; ðn þ m þ ‘ þ kÞ! k¼1
rk 2 D, (15a)
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
where mnk1 ¼
Z
1
W ðtÞtk1 pn ðtÞdt
(15b)
1
are the moments of pn ðtÞ. Taking into consideration the estimation (15), it is easily found that E nþmþ‘ ðjÞ ¼ 0 for jðtÞ 2 Pnþmþ‘1 , i.e. the integration formula is exact for polynomials of degree pn þ m þ ‘ 1. The Quadrature formula (11) without preassigned nodes, can also be obtained by extracting the singularities at fyi g‘i¼1 . In fact, we can write Z 1 ‘ Z 1 X wðtÞjðtÞ wðtÞjðtÞ dt ¼ dt, IðjÞ ¼ 0 ðy Þðt y Þ o ¼1 oðtÞ 1 i i i¼1 extracting the singularities, we have Z ‘ Z 1 ‘ X X wðtÞ½jðtÞ jðyi Þ jðyi Þ 1 wðtÞ dt þ IðjÞ ¼ dt. o0 ðyi Þðt yi Þ o0 ðyi Þ 1 t yi i¼1 1 i¼1 Using standard quadrature formula for the regular integral, it is obtained that: QpG nþ‘ ðjÞ
¼
‘ X n X i¼1 j¼1
‘ X jðtj Þ jðyi Þ jðyi Þ 2 c ðy Þ lj 0 0 ðy Þ 0 i o o ðyi Þðtj yi Þ i i¼1
(16a) and rearranging QpG nþ‘ ðjÞ
¼
n X j¼1
" # ‘ n X jðtj Þ X lj 2c0 ðyi Þ þ 0 lj jðyi Þ. o ðyi Þ oðtj Þ i¼1 j¼1 oðtj Þ (16b)
So
! n X lj 2c0 ðyi Þ n þ 0 li ¼ . o ðyi Þ oðtj Þ j¼1
761
1=Pnþm ðyi Þ which is of the order of 1=ri with ri ¼ tj yi (tj being the integration point nearest to yi ). This term takes its maximum value, which is 1=ii if tj coincides with xi ðyi ¼ xi þ ii Þ. In that case if i is very small, the error of the classical method becomes high. This could happen in the particular case where xi ¼ 0 and the number of integration points is odd ð2m þ 1; m ¼ 0; 1; 2; . . .Þ. Obviously in this case the error is maximum whereas in the case of even points the error decreases. The above qualitative estimation is obvious if we use analytical expressions for li , see for example the following relations (25), (27), (29), (33), (34) and (35). The estimation of the behaviour of such integrals can also be seen from the proposed test functions (36)–(39). The most important point is that, as n (the number of integration points) grows the situation worsens, i.e. as ri becomes smaller and smaller, cn ðyi Þ becomes higher and higher ðcn ðyi Þ behaves like yni Þ and consequently the error becomes very important. The conclusion is that if we use classical quadrature formulas, it is not clever to increase the number of integration points as the common practice is. 4. Some special cases of integrals with nearby poles In this paragraph we will give some more specific cases. Thus (i) we suppose that there are no preassigned nodes, i.e. OðtÞ ¼ 1. (ii) The weight is of the form, wðtÞ ¼ ð1 tÞa ð1 þ tÞb . (iii) The integration points are the roots of the corresponding orthogonal polynomial i.e. the Jacobi polynomial Pða;bÞ ðtÞ, and n Z 1 ða;bÞ ð1 tÞa ð1 þ tÞb Pða;bÞ ðtÞPða;bÞ dnm , n m ðtÞ dt ¼ hn 1
(16c)
The last formula can be also obtained directly from relation (13) again by extracting the singularities. If pn ðtÞ are selected to be the orthogonal polynomials with respect W(t) then the moments mk are equal to zero for kpn 1, and consequently taking into consideration (15),
21k Gðn þ a þ 1ÞGðn þ b þ 1Þ , ð2n k þ 1Þn!Gðn k þ 1Þ k ¼ ða þ bÞ; a; b4 1
hða;bÞ n
¼
ð17Þ Pða;bÞ ðtÞ n
or the roots of the orthogonal polynomial respect to the weight wðtÞ=oðtÞ. Z 1 ð1 tÞa ð1 þ tÞb ða;bÞ k Pn ðtÞt dt ¼ 0; kpn 1. oðtÞ 1
E nþmþ‘ ðjÞ ¼ 0 for jðtÞ 2 P2nþmþ‘1 .
4.1. Regular integrals with a single pole
Formula (11) is not a Gauss formula but it can be considered as a ‘‘modified or pseudo-Gauss’’ formula or as a Gauss formula with correction terms. A true Gauss formula must be exact for polynomials of degree p2n þ m þ 2‘ 1, since ðn þ m þ ‘Þ integration points are retained. The correction term (the second sum) in relation (11) represents the error created by the nearby singularities yi when a classical quadrature formula is used for the calculation of such integrals. The correction term, since cn ðyi Þ behaves like yni , is dominated by the behaviour of
The integral to be considered is Z 1 ð1 tÞa ð1 þ tÞb jðtÞ dt; Iðj; z0 Þ ¼ t z0 1
z0 eð1; 1Þ.
with
(18)
In this way the ‘‘modified or pseudo-Gauss’’ formula (11)–(13), takes the form QpG nþ1 ¼
n X lj jðtj Þ j¼1
z0 atj ,
t j z0
þ ln1 jðz0 Þ;
Pða;bÞ ðtj Þ ¼ 0, n ð19Þ
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
762
where lj is the classical Gauss weight given by lj ¼
2cða;bÞ ðtj Þ n
auxiliary results: ðicÞ ¼ cða;bÞ n
,
P0 ða;bÞ ðtj Þ n Z 1 1 ð1 tÞa ð1 þ tÞb Pða;bÞ ðtÞ n dt cða;bÞ ðzÞ ¼ n 2 1 zt
ð20Þ
and ln1 is the weight of the correction term which can be written as ! n X lj 2cða;bÞ ðz0 Þ ða;bÞ n n ¼ þ 2c0 ðz0 Þ . (21) l1 ¼ ða;bÞ t z0 Pn ðz0 Þ j¼1 j
4.2. Regular integrals with conjugate poles We consider the integral Z I 1 ðj; cÞ ¼
1
wðtÞjðtÞ dt ¼ 2 2 1 t þ c
Z
1
wðtÞjðtÞ dt ðt þ icÞðt icÞ 1
(22)
with poles at the points ic, where c51. Taking into consideration (6), (11) or (16), and considering the poles ic, as special integration points, it is finally obtained that: G CG QpG 1 ðj; cÞ ¼ Q1 ðj; cÞ þ Q1 ðj; cÞ
n X j¼1
QCG 1 ðj; cÞ ¼
lj
ðicÞ cða;bÞ n
Z ¼
n X lj Pn ðtÞ dt ¼ Pn ðicÞ t ic t ic j 1 j¼1 ! ic þ 1 , þ‘n ic 1 1
n X lj Pn ðtÞ dt ¼ Pn ðicÞ t þ ic t þ ic 1 j¼1 j ! ic 1 , þ‘n ic þ 1 1
ð24Þ
where lj and tj ðj ¼ 1; 2; . . . ; nÞ being the weights and the integration points respectively, of the classical GaussLegendre formula [47] arising from the classical GaussJacobi quadrature formula ða ¼ b ¼ 0Þ. Using (24), we obtain from (23) ( ) n n X X jðtj Þ lj 2 pG 1 1 Q1L ðj; cÞ ¼ tan lj 2 þ c c j¼1 t2j þ c2 tj þ c2 j¼1 ReðjðicÞÞ.
ð25Þ
The quadrature sum (25) is exact if jðtÞ is a polynomial of degree up to 2n þ 1.
(23a) 4.2.2. Chebyshev polynomials of (w(t) ¼ (1t2)1/2) If the weight w(t) is of the form
with QG 1 ðj; cÞ ¼
Z
jðtj Þ ; t2j þ c2
the
first
kind
wðtÞ ¼ ð1 t2 Þ1=2 ,
cða;bÞ ðicÞ cða;bÞ ðicÞ n n jðicÞ þ jðicÞ, ða;bÞ ða;bÞ icPn ðicÞ icPn ðicÞ ð23bÞ
where lj are determined from relation (20), ftj gnj¼1 are the zeros of the Jacobi polynomial Pða;bÞ ðtÞ corresponding to n w(t) [49] and cða;bÞ ðzÞ is also given by relation (20). The n degree of the polynomial accuracy of formula (23) is up to 2n+1. In (23) QG 1 constitutes the classical Gauss quadrature whereas QCG 1 , the additional term which corrects the influence of the nearby poles ic. As c become smaller the CG error in QG 1 becomes high and the significance of Q1 becomes important.
we denote by Tn(t) and Un1(t) the Chebyshev Polynomial of the first and second kind, respectively [45,46], which are orthogonal with respect the weights w(t) ¼ (1t2)1/2 and w(t) ¼ (1t2)1/2, respectively. The corresponding associated functions cða;bÞ ðzÞ, are defined [50], by n Z 2 ð1=2;1=2Þ 1 1 ð1 t2 Þ1=2 T n ðtÞ cn dt ðzÞ ¼ G n1 ðzÞ ¼ p p 1 tz 1 ¼ R0n ðzÞ ¼ ðz2 1Þ1=2 Rn ðzÞ, n na0; zeð1; 1Þ, Z 2 ð1=2;1=2Þ 1 1 ð1 t2 Þ1=2 U n1 ðtÞ cn dt ðzÞ ¼ Rn ðzÞ ¼ p p 1 zt ¼ ½z ðz2 1Þ1=2 n ; n40;
4.2.1. Legendre polynomials (w(t) ¼ 1) Relation (23) involves values of associated function cða;bÞ ðzÞ at ic. Taking into consideration relation (20), the n values of cða;bÞ ðzÞ at 7ic can be calculated exactly using the n quadrature formula (16) with special nodes at ic and ic, respectively. The quadrature formula (16) gives exact results because the integrand Pn(t) is a Legendre polynomial of degree n. Precisely, we have the following
zeð1; 1Þ. ð26Þ
The above functions are holomorphic in the complex z-plane ðz ¼ x þ iyÞ, cut along y ¼ 0, jxjp1. On this cut, it is true that: 2 1=2 R U n1 ðxÞ ¼ einW , n ðxÞ ¼ T n ðxÞ ið1 x Þ x ¼ cos y, 2 1=2 Rn ðxÞ. G n1 ðxÞ ¼ ið1 x Þ
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
where lj and tj ðj ¼ 1; 2; . . . ; nÞ are the weights determined from (20) and the corresponding integration points ðPða;bÞ ðtj Þ ¼ 0Þ [49] respectively, w1 is the weight function n of the additional term depended on the modified quadrature rule and cða;bÞ ðxÞ is also given by (20). n
Thus T n ðxÞ ¼ ½Rþ n ðxÞ þ Rn ðxÞ ¼ cos nW, U n1 ðxÞ ¼ ið1 x2 Þ1=2 ½Rþ n ðxÞ Rn ðxÞ=2 ¼ ½G þ n ðxÞ þ G n ðxÞ=2 ¼ sin nW= sin W.
Then (23) takes the form: QpG 1Cf ðj; cÞ ¼
n X j¼1
jðtj Þ lj 2 þ t j þ c2
(
n X p 1 pffiffiffiffiffiffiffiffiffiffiffiffiffi lj 2 2 t j þ c2 c 1þc j¼1
ReðjðicÞÞ,
)
ð27Þ
4.3.1. Legendre polynomials (w(t) ¼ 1) Eq. (32) takes the form QpG 2L ðj; cÞ ¼
with
2j 1 p ; tj ¼ cos 2n
lj ¼
p ; n
j ¼ 1; 2; . . . ; n.
4.2.3. Chebyshev polynomials of the second (w(t) ¼ (1t2)1/2) Similarly from (23), we obtain that: ( rffiffiffiffiffiffiffiffiffiffiffiffiffi ! n X jðtj Þ 1 pG lj 2 þ p 1þ 21 Q1Cs ðj; cÞ ¼ c t j þ c2 j¼1 ) n X lj ReðjðicÞÞ, 2 t þ c2 j¼1 j
(28) kind
p jp 2 lj ¼ sin ; nþ1 nþ1
ð29Þ
j ¼ 1; 2; . . . ; n. (30)
The quadrature sums (27) and (29) are exact if jðtÞ is polynomial of degree up to 2n þ 1. 4.3. Cauchy singular integrals with conjugate poles The Cauchy type singular integral of the form I 2 ðj; cÞ ¼
wðtÞ
jðtÞ dt; ðt xÞðt2 þ c2 Þ
1oxo1,
(31)
with poles at the points ic, ðc51Þ. Taking into consideration that ‘ X 1 1 1 ¼ þ , ðt xÞoðtÞ oðxÞðt xÞ i¼1 o0 ðyi Þðyi xÞðt yi Þ
a Cauchy type singular integral and two Riemann-type integrals, calculated according to the modified quadrature rule (11) which is proposed in the previous section, are derived. Finally, we get after some algebraic manipulations QpG 2 ðj; cÞ ¼
n X j¼1
(
lj
n X
jðtj Þ ðtj xÞðt2j þ c2 Þ j¼1 ( 1 Q ðxÞ jðxÞ 2 2 n 2 x þc Pn ðxÞ " # ) n X lj 2 1 1 þ tan Re½ðx þ icÞjðicÞ , c c j¼1 t2j þ c2 lj
ð33Þ
with jp ; tj ¼ cos nþ1
763
jðtj Þ 1 ðtj xÞðt2j þ c2 Þ x2 þ c2
) cða;bÞ ðxÞ n 2 ða;bÞ jðxÞ þ w1 Re½ðx þ icÞjðicÞ , Pn ðxÞ ð32Þ
where lj and tj, ðj ¼ 1; 2; . . . ; nÞ are the weights and the integration points of the Gauss–Legendre integration rule arising from the classical Gauss–Jacobi quadrature formula ða ¼ b ¼ 0Þ, Qn is the Legendre function of second kind [46,47] and Pn is the Legendre polynomial of degree n. The quadrature sum (33) is exact for polynomials of degree up to 2n þ 2. 4.3.2. Chebyshev polynomials of the first kind (w(t) ¼ (1t2)1/2) Taking into consideration (26), Eq. (32), takes the form n X jðtj Þ 1 H lj þ Q2Cf ðj; cÞ ¼ x2 þ c2 ðtj xÞðt2j þ c2 Þ j¼1 " ( U n1 ðxÞ p pffiffiffiffiffiffiffiffiffiffiffiffiffi jðxÞ p T n ðxÞ c 1 þ c2 ) # n X lj Re½ðx þ icÞjðicÞ , ð34Þ t2 þ c2 j¼1 j where lj and tj are determined from (28). 4.3.3. Chebyshev polynomial of the second kind (w(t)) ¼ (1t2)1/2) Taking into consideration (26), Eq. (32) takes the form QH 2Cs ðj; cÞ ¼
n X
jðtj Þ 1 2 2 þ c2 Þ x þ c2 xÞðt ðt j j j¼1 " ( rffiffiffiffiffiffiffiffiffiffiffiffiffi ! T nþ1 ðxÞ 1 jðxÞ þ p 1þ 21 p U n ðxÞ c ) # n X lj Re½ðx þ icÞjðicÞ , ð35Þ t 2 þ c2 j¼1 j lj
where lj and tj are determined by (30). The quadrature sums (34) and (35) are exact for polynomials of degree up to 2n þ 2.
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
764
5. Numerical applications The proposed quadrature formulas have been applied to four test functions and to the solution of a singular integral equation of the second kind, in order to have the opportunity to compare the modified rule with the ‘‘classical’’ Gauss formula [49], in the case that c ¼ 0:01. Test function (i), is the evaluation of the following integral: pffiffiffiffiffiffiffiffiffiffiffiffiffi Z 1 pffiffiffiffiffiffiffiffiffiffiffiffi2 1t c2 þ 1 c I 1Cs ¼ ffi 311:033380. dt ¼ p 2 2 c 1 t þ c (36) The integration points and the weights of the modified Gauss Chebyshev quadrature formula of the second kind (29), are determined from relations (30). From Table 1, it is observed that the modified Gauss–Chebyshev quadrature formula gives very good results for n ¼ 2. Test function (ii), is the evaluation of the following Cauchy-type singular integral: dt ðt xÞðt2 þ c2 Þ ‘njð1 xÞ=ð1 þ xÞj 2xðtan1 ð1=cÞÞ=c . ¼ x2 þ c2
I 2L ¼
I 2Cf ¼ 20941:15887 þ 20313:12297 ¼ 628:03590 that is the exact result for (38). It is observed that the deviation from the exact results, in the classical rule, may be very large. Test function (iv), is the evaluation of the following Cauchy-type singular integral: rffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi 1 px 1 2 dt ¼ 2 I 2Cs ¼ 1 þ 2. 1t 2 2 2 ðt xÞðt þ c Þ x þc c (39) The ‘‘classical’’ and the modified Gauss–Chebyshev quadrature rule of the second kind (relation (35)), for n ¼ 2 and x ¼ 2, give the exact result I 2Cs ¼ 0, whereas for n¼3 and x ¼ cos p=8, the ‘‘classical’’ gives, I 2Cs ffi 8505:19233, and the modified, I 2Cs ffi 8505:19233 þ 8165:171628 ¼ 340:0207,
ð37Þ
The ‘‘classical’’ and the modified Cauchy-type, Gauss– Legendre quadrature formula (relation (33)) for n ¼ 3 and x ¼ 0 give the exact value I 2L ¼ 0, whereas for n ¼ 3 and x ¼ :967802ðQn ðxÞ ¼ 0Þ, the classical gives, I 2L ¼ 9189:936960, and the modified: I 2L ffi 9189:93696 þ 8863:036297 ¼ 326:90663 that is the exact result. Hence with the modified rule we may have the correct value of (37). We may have the exact value of (37) and for n ¼ 1, because the exactness of the modified rule is up to 2n þ 2. Test function (iii), is the evaluation of the following Cauchy-type singular integral: I 2Cf
The ‘‘classical’’ and the modified Gauss–Chebyshev quadrature rules of the first kind (relation (34)), for n ¼ 2 and x ¼ 0, give the exact result I 2Cf ¼ 0, whereas for n ¼ 3 and x ¼ cos p=3ð¼ 0:5Þ, the ‘‘classical’’ gives, I 2Cf ¼ 20941:15887, and the modified
which is the exact value of (39). Finally, consider a singular integral equation of the second kind of the form Z 1 oðtÞ l dt þ kðt; xÞoðtÞdt ¼ px 2 þ1l ; tx x þ c2 1 c l ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi ð40Þ 1 þ c2 supplemented by the additional condition Z 1 oðtÞ dt ¼ pð1 lÞ,
(41)
1
where oðtÞ ¼ wðtÞj1 ðtÞ;
wðtÞ ¼ ð1 t2 Þ1=2 ;
kðt; xÞ ¼ t þ x.
Following Muskhelishvili [15] or Tsamasphyros and Theocaris [33], the expected solution of (40) and (41) is of the form.
1 1 pffiffiffiffiffiffiffiffiffiffiffiffi dt ¼ 2 1 t ðt xÞðt2 þ c2 Þ xp pffiffiffiffiffiffiffiffiffiffiffiffiffi . ¼ 2 cðx þ c2 Þ 1 þ c2
ð38Þ
Table 1 Differences between the classical Gauss–Chebyshev quadrature formula of the second kind, and the modified formula (29) I 1Cs ðj; 0:01Þ ffi 311:033380 n
Classical Gauss–Chebyshev quadrature rule of the second kind
Modified Gauss–Chebyshev quadrature rule of the second kind
2 3 4
6.280673 7855.552116 695.353
311.033310 311.033380 311.033380
j1 ðtÞ ¼ j2 ðtÞ=ðt2 þ c2 Þ.
(42)
Equations (40) and (41) could represent the solution of an antiplane crack going to the interface and loaded by two opposite concentrated forces applied to the points ð0 þ icÞ and ð0 icÞ. For the numerical evaluation of (40) and (41) the modified Gauss–Chebyshev quadratures of the first kind have been used. Especially for the nearby ðc ¼ 0:01Þ Cauchy singular integral of (40), quadrature (34) has been used whereas for the regular ones of (40) and (41) quadrature (27). For n ¼ 3 the modified quadratures give j2 ðt ¼ cos p=6Þ ¼ 0:75 which is the exact solution of (40) and (41), whereas the ‘‘classical’’ Gauss–Chebyshev of the second kind for n ¼ 3 gives j2 ðt ¼ cos p=6Þ ¼ 15:25 and for n ¼ 10, j2 ðt ¼ cos p=6Þ ¼ 203:30.
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
6. Conclusions and discussion In our study a ‘‘modified’’ or ‘‘pseudo-Gauss’’ quadrature rule has been proposed. Based on this rule, quadratures formulas result, which ensure the exact calculation of integrals of the form Z 1 wðtÞjðtÞ dt z0 eð1; 1Þ, 1 t z0 Z 1 wðtÞjðtÞ dt; c51 2 2 1 t þ c and wðtÞ
jðtÞ dt ðt xÞðt2 þ c2 Þ
765
crack in an infinite matrix loaded by two concentrated forces at two conjugate points of the infinite body. Although the proposed examples contain only first or second-degree polynomials, the proposed methods can be applied to polynomials of any degree as it is shown from the numerical solution of the system of singular integral Eqs. (40) and (41) where the exact result was obtained using only 3 integration points. From the proposed numerical applications it follows that for the ‘‘classical’’ Gauss formulas is rather impossible to provide the exact results for integrals with nearby singularities, and thus the necessity of using modified quadrature formulas in the case where the poles are very close to the integration interval, is imposed.
c51; 1oxo1,
where jðtÞ is a regular function. The ‘‘modified’’ or ‘‘pseudo-Gauss’’ rule is derived from a Lagrangian interpolatory procedure, where integration points are also taken at the poles 7ic ðc51Þ, of the function. The accuracy of the proposed quadrature formulas is two times greater than any other quadrature formulas. The disadvantage of these formulas is that, they cannot be applied directly for the solution of integral equations, because of the introduction of extra unknowns at the poles 7ic and consequently of the introduction of two more equations at 7ic in the final system. The advantages of the proposed formulas (11)–(13), compared with other methods and procedures [32–35], is that a great polynomial accuracy may be achieved and that the integration points and the weights of the classical Gauss rule, may also be used in the modified formulas (11)–(13). The proposed modified formulas are very effective in the case that the regular integrals and the Cauchy singular integrals, are given by formulas (18), (22) and (31), respectively, where jðtÞ is a holomorphic function in the complex z-plane or in the worst case a function with removable singularities at 71. The modified formulas have been derived from a procedure based on preassigned nodes, that is similar to the technique for the production of Lobatto or Radau formulas. The proposed quadrature formulas introduced additional terms that correct the influence of the nearby poles. In case that the poles move away from the interval of orthogonality, either the additional terms or the modified weights become insignificant and the ‘‘classical’’ Gauss formulas may be applied. The proposed quadrature formulas may also be used in three-dimensional problems in the case that the surface integrals can be evaluated by repeated one-dimensional integrations. The derived quadrature formulas may be easily implemented into a standard computer code. The proposed, in Section 5, basic functions, may appear in problems of fracture mechanics, when the stress parameters are calculated near a boundary. In particular singular integral (39) multiplied by (1t2)1/2 related to the solution of a
References [1] Theocaris PS, Tsamasphyros G. On the solution of boundary-value problems in the plane elasticity for multiply-connected regions (I: The second boundary-value problem). Lett Appl Eng Sci 1975;3:167–76. [2] Fromme J, Golberg MA. Numerical solution of a class of integral equations arising in two-dimensional aerodynamics. In: Golberg MA, editor. Solution methods for integral equations, theory and applications. New York: Plenum; 1978. [3] Moss WF. The two-dimensional oscillating airfoil: a new implementation of the Galerking method. SIAM J Numer Anal 1983;20:391–9. [4] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniquestheory and applications in engineering. Berlin, Heidelberg, New York: Springer; 1984. [5] Ioakimidis NI. Application of finite-part integrals to the singular integral equations of crack problems in plane and three-dimensional elasticity. Acta Mech 1982;45:31–47. [6] Tsamasphyros G. Numerical methods for fracture parameters calculation. In: Carpinteri A, editor. Handbook of fatigue crack propagation in metallic structures. Amsterdam: Elsevier Science B.V.; 1994. [7] Theotokoglou EE, Tsamasphyros G. Exact expressions for twodimensional singular and finite part integrals applicable in solid mechanics. Eng Anal Bound Elem 1998;22:125–32. [8] Mukherjee S. CPV and HFP integrals and their applications in the boundary element method. Int J Solids Struct 2000;37:6623–34. [9] Cruse TA, Aithal R. Non-singular boundary integral equation implementation. Int J Numer Methods Eng 1993;36:237–54. [10] Krishnasamy G, Rizzo FJ, Liu YJ. Boundary integral equations for thin bodies. Int J Numer Methods Eng 1994;37:107–21. [11] Tsamasphyros G, Theocaris PS. Integral-equation solution for half planes bonded together or in contact and containing internal cracks or holes. Ing Arch 1983;53:225–41. [12] Theotokoglou EN, Tsamasphyros G. Integral equations for any configuration of curved cracks and holes in an elastic strip. Ing Arch 1987;57:3–15. [13] Dirgantara T, Aliabadi MH. Crack growth analysis of plates loaded by bending and tension using dual boundary element method. Int J Fract 2000;105:27–74. [14] Aliabadi MH, Martin D. Boundary element hyper-singular formulation for elasto-plastic contact problems. Int J Numer Methods Eng 2000;48:995–1014. [15] Muskhelishvili NI. Singular integral equations. Groningen, The Netherlands: P. Noordhoff Ltd.; 1953. [16] Sladek N, Sladek J, Tanaka M. Regularization of hypersingular and nearly singular integrals in the potential theory and elasticity. Int J Numer Methods Eng 1993;36:1609–28. [17] Liu YJ, Rudolphi TJ. New identities for fundamental solutions and their applications to non-singular boundary element formulations. Comput Mech 1999;24:286–92.
ARTICLE IN PRESS 766
E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
[18] Mukherjee S. Finite parts of singular and hypersingular integrals with irregular boundary source points. Eng Anal Bound Elem 2000;24: 767–76. [19] Mukherjee S, Chati MK, Shi XL. Evaluation of nearly singular integrals in boundary element contour and node methods for threedimensional linear elasticity. Int J Solids Struct 2000;37:7633–54. [20] Telles JCF. A self-adaptive co-ordinate transformation for efficient numerical evaluation of general boundary element integrals. Int J Numer Methods Eng 1987;24:959–73. [21] Krishnasamy G, Schmerr LW, Rudolphi TJ, Rizzo FJ. Hypersingular boundary integral equations: some applications in acoustic and elastic waves scattering. ASME J Appl Mech 1990;57:404–14. [22] Dumont NA. On the efficient numerical evaluation of integrals with complex singularity poles. Eng Anal Bound Elem 1994;13:155–68. [23] Yun BI. A composite transformation for numerical integration of singular integrals in the BEM. Int J Numer Methods Eng 2003;57: 1883–98. [24] Kim P, Choi UJ. Two trigonometric quadrature formulas for evaluating hypersingular integrals. Int J Numer Methods Eng 2003;56:469–86. [25] Huang Q, Cruse TA. Some notes on singular integral techniques in Boundary Element Analysis. Int J Numer Methods Eng 1993;36: 2643–51. [26] Johnston PR. Application of sigmoidal transformations to weakly singular and nearly-singular boundary element integrals. Int J Numer Methods Eng 1999;45:1333–48. [27] Johnston PR, Elliott D. Error estimation of quadrature rules for evaluating singular integrals in boundary element problems. Int J Numer Methods Eng 2000;48:949–62. [28] Ma H, Kamiya N. A general algorithm for the numerical evaluation of nearly singular boundary integrals of various orders for two- and three-dimensional elasticity. Comput Mech 2002;29:277–88. [29] Korneichuk AA. Quadrature formulae for singular integrals. Z Vychisl Mat I Mat Fiz 1964;4/4(suppl.):64–74 [in Russian]. [30] Paget DE, Elliott D. An algorithm for the numerical evaluation of certain Cauchy principal value integrals. Numer Math 1972;19: 373–85. [31] Hunter DB. Some Gauss-type formulae for the evaluation of Cauchy principal values of integrals. Numer Math 1972;19:419–24. [32] Chawla MM, Ramakrishnan TR. Modified Gauss–Jacobi quadrature formulas for the numerical evaluation of Cauchy type singular integrals. BIT 1974;14:14–21. [33] Tsamasphyros G, Theocaris PS. Numerical solution of systems of singular integral equations with variable coefficients. Report.
[34]
[35]
[36] [37]
[38]
[39]
[40] [41]
[42] [43] [44] [45] [46] [47] [48] [49] [50]
National Technical University of Athens, November 1976, also. Appl Anal 1979;9:37–52. Tsamasphyros G, Theocaris PS. Sur une methode generale de quadrature pour des integrals du type Cauchy. Symp. Appl. Mathem., Salonica August 1976, also (1980). Rev Roum Sci Techn, Ser Mech Appl 1980;25:839–56. Ioakimidis NI, Theocaris PS. On the numerical evaluation of Cauchy principal value integrals. Symp. Of Appl. Mathem., Salonica August 1976, also (1977). Rev Roum Sci Techn, Ser Mech Appl 1977;22:803–18. Monegato G. On polynomials orthogonal with respect to particular variable-signed weight functions. ZAMP 1980;31:459–555. Ninis T. Factors influencing the solution of singular integral equations, Master Thesis (in Greek), National Technical University of Athens, 1980. Strouboulis T. Problems in the solution of singular integral equations, Master thesis (in Greek), National Technical University of Athens, 1980. Ioakimidis NI. A direct mehod for the construction of Gaussian quadrature rules for Cauchy type and finite part integrals. L’ Anal Numer Theorie Approx 1983;12:131–40. Monegato G. Quadrature formulas for functions with poles near the interval of integration. Math Comput 1986;47:137–58. Tsamasphyros G. A study of factors influencing the solution of singular integral equations. Eng Fract Mech 1986;24: 567–78. Tsamasphyros G, Dimou G. A Gauss quadrature formula for Cauchy type integrals. Comput Mech 1989;4:137–48. Tsamasphyros G, Dimou G. Gauss quadrature rules for finite part integrals. Int J Numer Methods Eng 1990;30:13–26. Laurie DP. Computation of Gauss-type quadrature formulas. J Comput Appl Math 2001;127:201–17. Isaacson E, Keller HB. Analysis of numerical methods. London: Wiley; 1966. Hildebrand FB. Introduction to numerical analysis. 2nd ed. New York: McGraw-Hill; 1974. Szego G. Orthogonal polynomials. New York: American Mathematical Society; 1959. Guelfond AO. Calcul des differences finies. Paris: Dunod; 1963. Stroud A, Secret D. Gaussian quadrature formulas. Englewood Cliffs, N.J: Prentice-Hall; 1966. Gladwell GML, England AH. Orthogonal polynomial solutions to some mixed boundary value problems in elasticity theory. Q J Mech Appl Math 1977;30:175–85.
Engineering Analysis with Boundary Elements 30 (2006) 758–766 www.elsevier.com/locate/enganabound
A modified Gauss quadrature formula with special integration points for evaluation of Quasi-singular integrals E.E. Theotokoglou, G. Tsamasphyros Faculty of Applied Sciences, Department of Mechanics, Laboratory of Strength of Materials, The National Technical Universtiy of Athens, Zoghraphou Campus, Theocaris Bld., GR-0157 73, Athens, Greece Received 11 November 2005; accepted 13 May 2006 Available online 14 July 2006
Abstract It is well known in the boundary element method that integration rules fail when the integrand presents a nearby singularity. This drawback arises when the field point is near the source point, i.e. in the case of a domain with very narrow boundaries or when the field point where we try to calculate stresses or any other field variables, is near the boundaries. In the present paper a quadrature formulas for ‘ isolated singularities near the integration interval, based on ordinary or special Langrange interpolatory polynomials, is obtained. This interpolatory formulas present similarities with known formulas for the numerical evaluation of singular integrals. Quadrature formulas for regular and singular integrals with conjugate poles are also derived. Numerical examples are given and the proposed quadrature rules present the expected polynomial accuracy. r 2006 Elsevier Ltd. All rights reserved. Keywords: Nearby poles; Quadrature formula; Lagrangian interpolation; Cauchy type integrals; Elasticity; Boundary element method; Numerical integration
1. Introduction In the solution of two- and three-dimensional problems with the boundary element method (BEM) [1–15], the BEM is required to evaluate integrals where the integrand is of the form 0(‘n(r)), 0(1/r) or 0(1/r2) in which r is the distance between the source point and a field point. If the source point belongs to the integration interval we say that the integral has a (i) weak singularity if the kernel behaves like ‘nr or 1=rg with go1 and the integrals are defined in the usual way. (ii) A strong singularity if the kernels behave like 1=rg with gX1. Integrals of this type are diverging integrals and can be defined as Hadamard finite part integrals [4–8]. If g ¼ 2k þ 1, k ¼ 0; 1; 2; . . ., the corresponding integrals can be defined as Cauchy principal value integrals [6–8]. Corresponding author. Sofouli 33, Nea Smirni 17122, Athens, Greece. Tel.: +30 210 772 1303; fax: +30 210 7721302. E-mail address: [email protected] (E.E. Theotokoglou).
0955-7997/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2006.05.001
Up to now many methods are adopted for the calculation of such integrals. So the integral is calculated by splitting the boundary in smaller or bigger integration intervals. The numerical evaluation of such integrals is very crucial in the BEM, not only for singular integrals but also for the regular ones, in order to compute field variables and their derivatives. The kernels of integral forms present a similar behavior ð1=rg Þ giving primary or secondary variables (like displacements, strains and stresses in the case of elasticity). Taking into consideration that secondary variables are mth order derivatives of primary variable entered on the boundary integral equation, their behaviour is like 1=rg with gX1. If the source point does not belong to the integration interval the integral is defined as a regular one. But if the distance r, of the source point and the nearest field point is small, the integrand takes a great value. The difficulties arise from the fact that the integrands of integrals along the boundary vary rapidly with the distance when the field point approaches the source point. If r is very small this error is very important and surprisingly is much more important as the number of integration points increase. We
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
say that these integrals present ‘‘a nearby singularity’’ or a quasi-singularity. In engineering problems, these difficulties arise in the analysis of thin or shell-like bodies [9,10], in bodies with very narrow boundaries, i.e. in I-section beams or in the necking problem or in curved thin members, in cracks near boundaries [11–13], in contact problems [14], and in the case where the boundary integral is constructed in order to calculate displacements and stresses at internal points, near the bounding surface of the body. In the particular case of a formulation based on Muskhelishvili’s complex potentials [1,11,15] in a plane elastic body D bounded by two smooth non-intersecting arcs L1 and L2 ðL ¼ L1 [ L2 Þ, the complex stress vector rn ðz0 Þ at a source point z0, is given by rn ðz0 Þ ¼ snn ðz0 Þ þ isnt ðz0 Þ ¼ Fðz0 Þ þ Fðz0 Þ o d¯z0 n z0 F0 ðz0 Þ þ Cðz0 Þ ; z0 2 D, þ dz0
ð1Þ
where snn ðz0 Þ is the normal component of the stress vector along an arbitrary fictitious cut ðt tÞ, snt ðz0 Þ is the corresponding tangential stress, and Fðz0 Þ, Cðz0 Þ are the Muskhelishvili’s complex potentials. These complex potentials can be expressed in terms of a distribution of complex dislocations g(t) [1,6,11] along L1 and L2, by the following singular integrals Z Z 1 gðtÞ 1 gðtÞ Fðz0 Þ ¼ dt þ dt, (2) 2pi L1 t z0 2pi L2 t z0 Cðz0 Þ ¼
Z Z 1 gðtÞ 1 gðtÞ dt þ dt 2pi L1 t z0 2pi L2 t z0 Z Z ¯tgðtÞ ¯tgðtÞ þ dt þ dt. 2 2 L1 ðt z0 Þ L2 ðt z0 Þ
If z0 is very close to L1 or L2 these integrals behave like 1=r ¼ j1=ðt z0 Þj, i.e. the integrals have a nearby singularity. Let now L1 and L2 are very narrow. If z0 belongs to L1 the integrals in (2) along L2 behave like nearby singular integrals and vice versa. The general form of an integral with nearby and strong singularities can be put in the following form: Z 1 wðtÞjðtÞ dt, (3a) IðjÞ ¼ 1 oðtÞ with oðtÞ ¼
‘ Y ðt yi Þ;
yi ¼ xi þ ii ,
(3b)
759
length s of a part of the boundary or of the whole boundary of the body D, can be transformed to an integral along ð1; 1Þ. Among the numerical methods for evaluation of nearby singularities, regularization procedures [16–19] and procedures trying to avoid singularities [20–24] appear. On the other hand, procedures based on transformations, which remedy the singularities [25–28] and on procedure based on general quadrature rules [29–44] have also been proposed. Strong, weak and nearby singularities have already been confronted at early eighties by our team [34,37,38,41]. Especially nearby singularities are confronted in two master thesis [37,38] and the results are incorporated in our codes for the solution of boundary integral equations. In the present study, we have expanded the known Gauss-integration rules in order to confront the above problem where many strong or weak singularities can coexist. For the sake of generality we have considered the most general case by supposing that, for some reasons, we need to include some preassigned nodes. A ‘‘modified’’ or ‘‘pseudo-Gauss’’ quadrature formulas has been proposed. This quadrature formula is produced from a simple Lagrange interpolation jðtÞ. For question of generality preassigned nodes are considered, so Lobatto, Radau and Kronrod formulas [34,35,44] are derived. Based on the proposed rule, we have also constructed modified quadrature formulas for regular integrals with a single pole, for regular integrals with conjugate poles: Z 1 wðtÞjðtÞ dt; c51, 2 2 1 t þ c where the polynomial accuracy would be equal to 2n+1 for n integration points and for Cauchy singular integral with conjugate poles: Z 1 wðtÞjðtÞ dt; c51; 1oxo1, 2 2 1 ðt xÞðt þ c Þ where the polynomial accuracy would be equal to 2n+2. Thus, a greater polynomial accuracy has been achieved in the case of the proposed formula, in comparison with other methods and procedures [32–35]. Finally, examples are presented for the application of the modified quadrature rule illustrating, regular and singular Cauchy type integrals with conjugate poles. 2. Some definitions
i¼1
where ei is small or zero and xi is on ½1; 1 or very close to this interval and w(t) the weight function which is fixed positive and integrable on ½1; 1, given by wðtÞ ¼ ð1 tÞa ð1 þ tÞb ;
a; b4 1.
(3c)
The previous form (3) concerns not only line integrals but also surface integrals if the latter can be reduced to one-dimensional integrals and then the form (3) can be applied. Obviously any integral with respect to the arc
Let us give some necessary definitions: We denote by D a simply connected domain containing ð1; 1Þ in its interior. We denote by jðzÞ the analytic continuation of jðtÞ into D. Let C a simple contour inside D, presenting angular points at 1 and 1 with angles a and b, respectively (Fig. 1). These angular points are due to the behaviour of jðzÞ.
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
760
or in more compact form. jnþmþ‘ ðtÞ ¼
nþm X j¼1
þ
oðtÞPnþm ðtÞ jðtj Þ ðt tj Þoðtj ÞP0nþm ðtj Þ
‘ X i¼1
oðtÞPnþm ðtÞ jðyi Þ ðt yi Þo0 ðyi ÞPnþm ðyi Þ
ð8Þ
and Fig. 1. A simple connected domain D containing (1, 1), at its interior.
oðtÞPnþm ðtÞ Rnþmþ‘ ðtÞ ¼ 2pi
I
jðzÞ dz. ðz tÞoðzÞPnþm ðzÞ
(9)
C
ftj gnj¼1 ; fzk gm k¼1
Let two sets of points inside D, belonging or not to the integration interval. We suppose that zk are preassigned, whereas tj can be selected arbitrarily. These points will be used as interpolation points. Besides these interpolation points, it is possible to use as ‘‘special integration points’’ the roots yi ði ¼ 1; 2; . . . ; ‘Þ of oðtÞ ¼ 0. In this case the interpolation polynomial jnþmþ‘ ðzÞ is a complex polynomial of degree ðn þ m þ ‘ 1Þ, and we must have jðtj Þ ¼ jnþmþ‘ ðtj Þ; where ( tj ¼
tj ; zk ;
jðyi Þ ¼ jnþmþ‘ ðyi Þ,
(4)
1pjpn; j ¼ n þ kpn þ m:
We denote by pn ðtÞ ¼ Y nþm
n Y
ðt tj Þ;
OðtÞ ¼
j¼1
m Y
ðt zk Þ,
k¼1
ðtÞ ¼ pn ðtÞOðtÞ;
W ðtÞ ¼ wðtÞOðtÞ.
ð5Þ
3. Quadrature formulas with special integration or collocation points (modified or pseudo-Gauss formula with correction terms) As we mentioned above, we may use as special integration points the roots yi ði ¼ 1; 2; . . . ; ‘Þ of oðtÞ. Thus from the interpolation formula (4) and taking into consideration Lagrange procedure [45–47], it is obtained jðtÞ ¼ jnþmþ‘ ðtÞ þ Rnþmþ‘ ðtÞ, jnþmþ‘ ðtÞ ¼
n X j¼1
þ þ
(6)
oðtÞPnþm ðtÞ jðtj Þ ðt tj Þoðtj ÞOðtj Þp0n ðtj Þ
n X
oðtÞPnþm ðtÞ jðzk Þ ðt zk Þoðzk ÞO0 ðzk Þpn ðzk Þ k¼1
‘ X i¼1
oðtÞPnþm ðtÞ jðyi Þ, ðt yi Þo0 ðyi ÞPnþm ðyi Þ
Multiplying both sides of (6) by wðtÞ=oðtÞ and integrating on ½1; 1, we get Z 1 wðtÞ jðtÞdt ¼ QpG (10) nþmþ‘ ðjÞ þ E nþmþ‘ ðjÞ, oðtÞ 1 where QpG nþmþ‘ ðjÞ ¼
lj
j¼1
‘ jðtj Þ X þ ln jðyi Þ oðtj Þ i¼1 i
is the quadrature formula, Z 1 2c ðtj Þ 1 W ðtÞpn ðtÞ lj ¼ 0 dt ¼ 0 n Pnþm ðtj Þ 1 t tj Pnþm ðtj Þ
(11)
(12)
are the classical weights of the standard interpolatory formula for regular integrals, Z 1 1 W ðtÞpn ðtÞ n li ¼ 0 dt o ðyi ÞPnþm ðyi Þ 1 t yi 2cn ðyi Þ ð13Þ ¼ 0 o ðyi ÞPnþm ðyi Þ are the weights corresponding to the special integration points, and Z 1 I 1 jðzÞ W ðtÞpn ðtÞ dt dz E nþmþ‘ ðjÞ ¼ 2pi oðzÞPnþm ðzÞ 1 z t C I 1 jðzÞcn ðzÞ ¼ dz ð14Þ pi oðzÞPnþm ðzÞ C
is the error of the quadrature formula. The first term of (11) can also be obtained if the classical integration formula is applied to the integrand jðtÞ=oðtÞ. If, in other words, the singularities (the poles) of the integrand are ignored, the second term can be considered as a ‘‘correction term’’. The error term following the prescribed procedure and Guelfond [48], can be written as E nþmþ‘ ðjÞ ¼
ð7Þ
nþm X
1 X
mnk1 jðnþmþ‘þkÞ ðrk Þ; ðn þ m þ ‘ þ kÞ! k¼1
rk 2 D, (15a)
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
where mnk1 ¼
Z
1
W ðtÞtk1 pn ðtÞdt
(15b)
1
are the moments of pn ðtÞ. Taking into consideration the estimation (15), it is easily found that E nþmþ‘ ðjÞ ¼ 0 for jðtÞ 2 Pnþmþ‘1 , i.e. the integration formula is exact for polynomials of degree pn þ m þ ‘ 1. The Quadrature formula (11) without preassigned nodes, can also be obtained by extracting the singularities at fyi g‘i¼1 . In fact, we can write Z 1 ‘ Z 1 X wðtÞjðtÞ wðtÞjðtÞ dt ¼ dt, IðjÞ ¼ 0 ðy Þðt y Þ o ¼1 oðtÞ 1 i i i¼1 extracting the singularities, we have Z ‘ Z 1 ‘ X X wðtÞ½jðtÞ jðyi Þ jðyi Þ 1 wðtÞ dt þ IðjÞ ¼ dt. o0 ðyi Þðt yi Þ o0 ðyi Þ 1 t yi i¼1 1 i¼1 Using standard quadrature formula for the regular integral, it is obtained that: QpG nþ‘ ðjÞ
¼
‘ X n X i¼1 j¼1
‘ X jðtj Þ jðyi Þ jðyi Þ 2 c ðy Þ lj 0 0 ðy Þ 0 i o o ðyi Þðtj yi Þ i i¼1
(16a) and rearranging QpG nþ‘ ðjÞ
¼
n X j¼1
" # ‘ n X jðtj Þ X lj 2c0 ðyi Þ þ 0 lj jðyi Þ. o ðyi Þ oðtj Þ i¼1 j¼1 oðtj Þ (16b)
So
! n X lj 2c0 ðyi Þ n þ 0 li ¼ . o ðyi Þ oðtj Þ j¼1
761
1=Pnþm ðyi Þ which is of the order of 1=ri with ri ¼ tj yi (tj being the integration point nearest to yi ). This term takes its maximum value, which is 1=ii if tj coincides with xi ðyi ¼ xi þ ii Þ. In that case if i is very small, the error of the classical method becomes high. This could happen in the particular case where xi ¼ 0 and the number of integration points is odd ð2m þ 1; m ¼ 0; 1; 2; . . .Þ. Obviously in this case the error is maximum whereas in the case of even points the error decreases. The above qualitative estimation is obvious if we use analytical expressions for li , see for example the following relations (25), (27), (29), (33), (34) and (35). The estimation of the behaviour of such integrals can also be seen from the proposed test functions (36)–(39). The most important point is that, as n (the number of integration points) grows the situation worsens, i.e. as ri becomes smaller and smaller, cn ðyi Þ becomes higher and higher ðcn ðyi Þ behaves like yni Þ and consequently the error becomes very important. The conclusion is that if we use classical quadrature formulas, it is not clever to increase the number of integration points as the common practice is. 4. Some special cases of integrals with nearby poles In this paragraph we will give some more specific cases. Thus (i) we suppose that there are no preassigned nodes, i.e. OðtÞ ¼ 1. (ii) The weight is of the form, wðtÞ ¼ ð1 tÞa ð1 þ tÞb . (iii) The integration points are the roots of the corresponding orthogonal polynomial i.e. the Jacobi polynomial Pða;bÞ ðtÞ, and n Z 1 ða;bÞ ð1 tÞa ð1 þ tÞb Pða;bÞ ðtÞPða;bÞ dnm , n m ðtÞ dt ¼ hn 1
(16c)
The last formula can be also obtained directly from relation (13) again by extracting the singularities. If pn ðtÞ are selected to be the orthogonal polynomials with respect W(t) then the moments mk are equal to zero for kpn 1, and consequently taking into consideration (15),
21k Gðn þ a þ 1ÞGðn þ b þ 1Þ , ð2n k þ 1Þn!Gðn k þ 1Þ k ¼ ða þ bÞ; a; b4 1
hða;bÞ n
¼
ð17Þ Pða;bÞ ðtÞ n
or the roots of the orthogonal polynomial respect to the weight wðtÞ=oðtÞ. Z 1 ð1 tÞa ð1 þ tÞb ða;bÞ k Pn ðtÞt dt ¼ 0; kpn 1. oðtÞ 1
E nþmþ‘ ðjÞ ¼ 0 for jðtÞ 2 P2nþmþ‘1 .
4.1. Regular integrals with a single pole
Formula (11) is not a Gauss formula but it can be considered as a ‘‘modified or pseudo-Gauss’’ formula or as a Gauss formula with correction terms. A true Gauss formula must be exact for polynomials of degree p2n þ m þ 2‘ 1, since ðn þ m þ ‘Þ integration points are retained. The correction term (the second sum) in relation (11) represents the error created by the nearby singularities yi when a classical quadrature formula is used for the calculation of such integrals. The correction term, since cn ðyi Þ behaves like yni , is dominated by the behaviour of
The integral to be considered is Z 1 ð1 tÞa ð1 þ tÞb jðtÞ dt; Iðj; z0 Þ ¼ t z0 1
z0 eð1; 1Þ.
with
(18)
In this way the ‘‘modified or pseudo-Gauss’’ formula (11)–(13), takes the form QpG nþ1 ¼
n X lj jðtj Þ j¼1
z0 atj ,
t j z0
þ ln1 jðz0 Þ;
Pða;bÞ ðtj Þ ¼ 0, n ð19Þ
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
762
where lj is the classical Gauss weight given by lj ¼
2cða;bÞ ðtj Þ n
auxiliary results: ðicÞ ¼ cða;bÞ n
,
P0 ða;bÞ ðtj Þ n Z 1 1 ð1 tÞa ð1 þ tÞb Pða;bÞ ðtÞ n dt cða;bÞ ðzÞ ¼ n 2 1 zt
ð20Þ
and ln1 is the weight of the correction term which can be written as ! n X lj 2cða;bÞ ðz0 Þ ða;bÞ n n ¼ þ 2c0 ðz0 Þ . (21) l1 ¼ ða;bÞ t z0 Pn ðz0 Þ j¼1 j
4.2. Regular integrals with conjugate poles We consider the integral Z I 1 ðj; cÞ ¼
1
wðtÞjðtÞ dt ¼ 2 2 1 t þ c
Z
1
wðtÞjðtÞ dt ðt þ icÞðt icÞ 1
(22)
with poles at the points ic, where c51. Taking into consideration (6), (11) or (16), and considering the poles ic, as special integration points, it is finally obtained that: G CG QpG 1 ðj; cÞ ¼ Q1 ðj; cÞ þ Q1 ðj; cÞ
n X j¼1
QCG 1 ðj; cÞ ¼
lj
ðicÞ cða;bÞ n
Z ¼
n X lj Pn ðtÞ dt ¼ Pn ðicÞ t ic t ic j 1 j¼1 ! ic þ 1 , þ‘n ic 1 1
n X lj Pn ðtÞ dt ¼ Pn ðicÞ t þ ic t þ ic 1 j¼1 j ! ic 1 , þ‘n ic þ 1 1
ð24Þ
where lj and tj ðj ¼ 1; 2; . . . ; nÞ being the weights and the integration points respectively, of the classical GaussLegendre formula [47] arising from the classical GaussJacobi quadrature formula ða ¼ b ¼ 0Þ. Using (24), we obtain from (23) ( ) n n X X jðtj Þ lj 2 pG 1 1 Q1L ðj; cÞ ¼ tan lj 2 þ c c j¼1 t2j þ c2 tj þ c2 j¼1 ReðjðicÞÞ.
ð25Þ
The quadrature sum (25) is exact if jðtÞ is a polynomial of degree up to 2n þ 1.
(23a) 4.2.2. Chebyshev polynomials of (w(t) ¼ (1t2)1/2) If the weight w(t) is of the form
with QG 1 ðj; cÞ ¼
Z
jðtj Þ ; t2j þ c2
the
first
kind
wðtÞ ¼ ð1 t2 Þ1=2 ,
cða;bÞ ðicÞ cða;bÞ ðicÞ n n jðicÞ þ jðicÞ, ða;bÞ ða;bÞ icPn ðicÞ icPn ðicÞ ð23bÞ
where lj are determined from relation (20), ftj gnj¼1 are the zeros of the Jacobi polynomial Pða;bÞ ðtÞ corresponding to n w(t) [49] and cða;bÞ ðzÞ is also given by relation (20). The n degree of the polynomial accuracy of formula (23) is up to 2n+1. In (23) QG 1 constitutes the classical Gauss quadrature whereas QCG 1 , the additional term which corrects the influence of the nearby poles ic. As c become smaller the CG error in QG 1 becomes high and the significance of Q1 becomes important.
we denote by Tn(t) and Un1(t) the Chebyshev Polynomial of the first and second kind, respectively [45,46], which are orthogonal with respect the weights w(t) ¼ (1t2)1/2 and w(t) ¼ (1t2)1/2, respectively. The corresponding associated functions cða;bÞ ðzÞ, are defined [50], by n Z 2 ð1=2;1=2Þ 1 1 ð1 t2 Þ1=2 T n ðtÞ cn dt ðzÞ ¼ G n1 ðzÞ ¼ p p 1 tz 1 ¼ R0n ðzÞ ¼ ðz2 1Þ1=2 Rn ðzÞ, n na0; zeð1; 1Þ, Z 2 ð1=2;1=2Þ 1 1 ð1 t2 Þ1=2 U n1 ðtÞ cn dt ðzÞ ¼ Rn ðzÞ ¼ p p 1 zt ¼ ½z ðz2 1Þ1=2 n ; n40;
4.2.1. Legendre polynomials (w(t) ¼ 1) Relation (23) involves values of associated function cða;bÞ ðzÞ at ic. Taking into consideration relation (20), the n values of cða;bÞ ðzÞ at 7ic can be calculated exactly using the n quadrature formula (16) with special nodes at ic and ic, respectively. The quadrature formula (16) gives exact results because the integrand Pn(t) is a Legendre polynomial of degree n. Precisely, we have the following
zeð1; 1Þ. ð26Þ
The above functions are holomorphic in the complex z-plane ðz ¼ x þ iyÞ, cut along y ¼ 0, jxjp1. On this cut, it is true that: 2 1=2 R U n1 ðxÞ ¼ einW , n ðxÞ ¼ T n ðxÞ ið1 x Þ x ¼ cos y, 2 1=2 Rn ðxÞ. G n1 ðxÞ ¼ ið1 x Þ
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
where lj and tj ðj ¼ 1; 2; . . . ; nÞ are the weights determined from (20) and the corresponding integration points ðPða;bÞ ðtj Þ ¼ 0Þ [49] respectively, w1 is the weight function n of the additional term depended on the modified quadrature rule and cða;bÞ ðxÞ is also given by (20). n
Thus T n ðxÞ ¼ ½Rþ n ðxÞ þ Rn ðxÞ ¼ cos nW, U n1 ðxÞ ¼ ið1 x2 Þ1=2 ½Rþ n ðxÞ Rn ðxÞ=2 ¼ ½G þ n ðxÞ þ G n ðxÞ=2 ¼ sin nW= sin W.
Then (23) takes the form: QpG 1Cf ðj; cÞ ¼
n X j¼1
jðtj Þ lj 2 þ t j þ c2
(
n X p 1 pffiffiffiffiffiffiffiffiffiffiffiffiffi lj 2 2 t j þ c2 c 1þc j¼1
ReðjðicÞÞ,
)
ð27Þ
4.3.1. Legendre polynomials (w(t) ¼ 1) Eq. (32) takes the form QpG 2L ðj; cÞ ¼
with
2j 1 p ; tj ¼ cos 2n
lj ¼
p ; n
j ¼ 1; 2; . . . ; n.
4.2.3. Chebyshev polynomials of the second (w(t) ¼ (1t2)1/2) Similarly from (23), we obtain that: ( rffiffiffiffiffiffiffiffiffiffiffiffiffi ! n X jðtj Þ 1 pG lj 2 þ p 1þ 21 Q1Cs ðj; cÞ ¼ c t j þ c2 j¼1 ) n X lj ReðjðicÞÞ, 2 t þ c2 j¼1 j
(28) kind
p jp 2 lj ¼ sin ; nþ1 nþ1
ð29Þ
j ¼ 1; 2; . . . ; n. (30)
The quadrature sums (27) and (29) are exact if jðtÞ is polynomial of degree up to 2n þ 1. 4.3. Cauchy singular integrals with conjugate poles The Cauchy type singular integral of the form I 2 ðj; cÞ ¼
wðtÞ
jðtÞ dt; ðt xÞðt2 þ c2 Þ
1oxo1,
(31)
with poles at the points ic, ðc51Þ. Taking into consideration that ‘ X 1 1 1 ¼ þ , ðt xÞoðtÞ oðxÞðt xÞ i¼1 o0 ðyi Þðyi xÞðt yi Þ
a Cauchy type singular integral and two Riemann-type integrals, calculated according to the modified quadrature rule (11) which is proposed in the previous section, are derived. Finally, we get after some algebraic manipulations QpG 2 ðj; cÞ ¼
n X j¼1
(
lj
n X
jðtj Þ ðtj xÞðt2j þ c2 Þ j¼1 ( 1 Q ðxÞ jðxÞ 2 2 n 2 x þc Pn ðxÞ " # ) n X lj 2 1 1 þ tan Re½ðx þ icÞjðicÞ , c c j¼1 t2j þ c2 lj
ð33Þ
with jp ; tj ¼ cos nþ1
763
jðtj Þ 1 ðtj xÞðt2j þ c2 Þ x2 þ c2
) cða;bÞ ðxÞ n 2 ða;bÞ jðxÞ þ w1 Re½ðx þ icÞjðicÞ , Pn ðxÞ ð32Þ
where lj and tj, ðj ¼ 1; 2; . . . ; nÞ are the weights and the integration points of the Gauss–Legendre integration rule arising from the classical Gauss–Jacobi quadrature formula ða ¼ b ¼ 0Þ, Qn is the Legendre function of second kind [46,47] and Pn is the Legendre polynomial of degree n. The quadrature sum (33) is exact for polynomials of degree up to 2n þ 2. 4.3.2. Chebyshev polynomials of the first kind (w(t) ¼ (1t2)1/2) Taking into consideration (26), Eq. (32), takes the form n X jðtj Þ 1 H lj þ Q2Cf ðj; cÞ ¼ x2 þ c2 ðtj xÞðt2j þ c2 Þ j¼1 " ( U n1 ðxÞ p pffiffiffiffiffiffiffiffiffiffiffiffiffi jðxÞ p T n ðxÞ c 1 þ c2 ) # n X lj Re½ðx þ icÞjðicÞ , ð34Þ t2 þ c2 j¼1 j where lj and tj are determined from (28). 4.3.3. Chebyshev polynomial of the second kind (w(t)) ¼ (1t2)1/2) Taking into consideration (26), Eq. (32) takes the form QH 2Cs ðj; cÞ ¼
n X
jðtj Þ 1 2 2 þ c2 Þ x þ c2 xÞðt ðt j j j¼1 " ( rffiffiffiffiffiffiffiffiffiffiffiffiffi ! T nþ1 ðxÞ 1 jðxÞ þ p 1þ 21 p U n ðxÞ c ) # n X lj Re½ðx þ icÞjðicÞ , ð35Þ t 2 þ c2 j¼1 j lj
where lj and tj are determined by (30). The quadrature sums (34) and (35) are exact for polynomials of degree up to 2n þ 2.
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
764
5. Numerical applications The proposed quadrature formulas have been applied to four test functions and to the solution of a singular integral equation of the second kind, in order to have the opportunity to compare the modified rule with the ‘‘classical’’ Gauss formula [49], in the case that c ¼ 0:01. Test function (i), is the evaluation of the following integral: pffiffiffiffiffiffiffiffiffiffiffiffiffi Z 1 pffiffiffiffiffiffiffiffiffiffiffiffi2 1t c2 þ 1 c I 1Cs ¼ ffi 311:033380. dt ¼ p 2 2 c 1 t þ c (36) The integration points and the weights of the modified Gauss Chebyshev quadrature formula of the second kind (29), are determined from relations (30). From Table 1, it is observed that the modified Gauss–Chebyshev quadrature formula gives very good results for n ¼ 2. Test function (ii), is the evaluation of the following Cauchy-type singular integral: dt ðt xÞðt2 þ c2 Þ ‘njð1 xÞ=ð1 þ xÞj 2xðtan1 ð1=cÞÞ=c . ¼ x2 þ c2
I 2L ¼
I 2Cf ¼ 20941:15887 þ 20313:12297 ¼ 628:03590 that is the exact result for (38). It is observed that the deviation from the exact results, in the classical rule, may be very large. Test function (iv), is the evaluation of the following Cauchy-type singular integral: rffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi 1 px 1 2 dt ¼ 2 I 2Cs ¼ 1 þ 2. 1t 2 2 2 ðt xÞðt þ c Þ x þc c (39) The ‘‘classical’’ and the modified Gauss–Chebyshev quadrature rule of the second kind (relation (35)), for n ¼ 2 and x ¼ 2, give the exact result I 2Cs ¼ 0, whereas for n¼3 and x ¼ cos p=8, the ‘‘classical’’ gives, I 2Cs ffi 8505:19233, and the modified, I 2Cs ffi 8505:19233 þ 8165:171628 ¼ 340:0207,
ð37Þ
The ‘‘classical’’ and the modified Cauchy-type, Gauss– Legendre quadrature formula (relation (33)) for n ¼ 3 and x ¼ 0 give the exact value I 2L ¼ 0, whereas for n ¼ 3 and x ¼ :967802ðQn ðxÞ ¼ 0Þ, the classical gives, I 2L ¼ 9189:936960, and the modified: I 2L ffi 9189:93696 þ 8863:036297 ¼ 326:90663 that is the exact result. Hence with the modified rule we may have the correct value of (37). We may have the exact value of (37) and for n ¼ 1, because the exactness of the modified rule is up to 2n þ 2. Test function (iii), is the evaluation of the following Cauchy-type singular integral: I 2Cf
The ‘‘classical’’ and the modified Gauss–Chebyshev quadrature rules of the first kind (relation (34)), for n ¼ 2 and x ¼ 0, give the exact result I 2Cf ¼ 0, whereas for n ¼ 3 and x ¼ cos p=3ð¼ 0:5Þ, the ‘‘classical’’ gives, I 2Cf ¼ 20941:15887, and the modified
which is the exact value of (39). Finally, consider a singular integral equation of the second kind of the form Z 1 oðtÞ l dt þ kðt; xÞoðtÞdt ¼ px 2 þ1l ; tx x þ c2 1 c l ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi ð40Þ 1 þ c2 supplemented by the additional condition Z 1 oðtÞ dt ¼ pð1 lÞ,
(41)
1
where oðtÞ ¼ wðtÞj1 ðtÞ;
wðtÞ ¼ ð1 t2 Þ1=2 ;
kðt; xÞ ¼ t þ x.
Following Muskhelishvili [15] or Tsamasphyros and Theocaris [33], the expected solution of (40) and (41) is of the form.
1 1 pffiffiffiffiffiffiffiffiffiffiffiffi dt ¼ 2 1 t ðt xÞðt2 þ c2 Þ xp pffiffiffiffiffiffiffiffiffiffiffiffiffi . ¼ 2 cðx þ c2 Þ 1 þ c2
ð38Þ
Table 1 Differences between the classical Gauss–Chebyshev quadrature formula of the second kind, and the modified formula (29) I 1Cs ðj; 0:01Þ ffi 311:033380 n
Classical Gauss–Chebyshev quadrature rule of the second kind
Modified Gauss–Chebyshev quadrature rule of the second kind
2 3 4
6.280673 7855.552116 695.353
311.033310 311.033380 311.033380
j1 ðtÞ ¼ j2 ðtÞ=ðt2 þ c2 Þ.
(42)
Equations (40) and (41) could represent the solution of an antiplane crack going to the interface and loaded by two opposite concentrated forces applied to the points ð0 þ icÞ and ð0 icÞ. For the numerical evaluation of (40) and (41) the modified Gauss–Chebyshev quadratures of the first kind have been used. Especially for the nearby ðc ¼ 0:01Þ Cauchy singular integral of (40), quadrature (34) has been used whereas for the regular ones of (40) and (41) quadrature (27). For n ¼ 3 the modified quadratures give j2 ðt ¼ cos p=6Þ ¼ 0:75 which is the exact solution of (40) and (41), whereas the ‘‘classical’’ Gauss–Chebyshev of the second kind for n ¼ 3 gives j2 ðt ¼ cos p=6Þ ¼ 15:25 and for n ¼ 10, j2 ðt ¼ cos p=6Þ ¼ 203:30.
ARTICLE IN PRESS E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
6. Conclusions and discussion In our study a ‘‘modified’’ or ‘‘pseudo-Gauss’’ quadrature rule has been proposed. Based on this rule, quadratures formulas result, which ensure the exact calculation of integrals of the form Z 1 wðtÞjðtÞ dt z0 eð1; 1Þ, 1 t z0 Z 1 wðtÞjðtÞ dt; c51 2 2 1 t þ c and wðtÞ
jðtÞ dt ðt xÞðt2 þ c2 Þ
765
crack in an infinite matrix loaded by two concentrated forces at two conjugate points of the infinite body. Although the proposed examples contain only first or second-degree polynomials, the proposed methods can be applied to polynomials of any degree as it is shown from the numerical solution of the system of singular integral Eqs. (40) and (41) where the exact result was obtained using only 3 integration points. From the proposed numerical applications it follows that for the ‘‘classical’’ Gauss formulas is rather impossible to provide the exact results for integrals with nearby singularities, and thus the necessity of using modified quadrature formulas in the case where the poles are very close to the integration interval, is imposed.
c51; 1oxo1,
where jðtÞ is a regular function. The ‘‘modified’’ or ‘‘pseudo-Gauss’’ rule is derived from a Lagrangian interpolatory procedure, where integration points are also taken at the poles 7ic ðc51Þ, of the function. The accuracy of the proposed quadrature formulas is two times greater than any other quadrature formulas. The disadvantage of these formulas is that, they cannot be applied directly for the solution of integral equations, because of the introduction of extra unknowns at the poles 7ic and consequently of the introduction of two more equations at 7ic in the final system. The advantages of the proposed formulas (11)–(13), compared with other methods and procedures [32–35], is that a great polynomial accuracy may be achieved and that the integration points and the weights of the classical Gauss rule, may also be used in the modified formulas (11)–(13). The proposed modified formulas are very effective in the case that the regular integrals and the Cauchy singular integrals, are given by formulas (18), (22) and (31), respectively, where jðtÞ is a holomorphic function in the complex z-plane or in the worst case a function with removable singularities at 71. The modified formulas have been derived from a procedure based on preassigned nodes, that is similar to the technique for the production of Lobatto or Radau formulas. The proposed quadrature formulas introduced additional terms that correct the influence of the nearby poles. In case that the poles move away from the interval of orthogonality, either the additional terms or the modified weights become insignificant and the ‘‘classical’’ Gauss formulas may be applied. The proposed quadrature formulas may also be used in three-dimensional problems in the case that the surface integrals can be evaluated by repeated one-dimensional integrations. The derived quadrature formulas may be easily implemented into a standard computer code. The proposed, in Section 5, basic functions, may appear in problems of fracture mechanics, when the stress parameters are calculated near a boundary. In particular singular integral (39) multiplied by (1t2)1/2 related to the solution of a
References [1] Theocaris PS, Tsamasphyros G. On the solution of boundary-value problems in the plane elasticity for multiply-connected regions (I: The second boundary-value problem). Lett Appl Eng Sci 1975;3:167–76. [2] Fromme J, Golberg MA. Numerical solution of a class of integral equations arising in two-dimensional aerodynamics. In: Golberg MA, editor. Solution methods for integral equations, theory and applications. New York: Plenum; 1978. [3] Moss WF. The two-dimensional oscillating airfoil: a new implementation of the Galerking method. SIAM J Numer Anal 1983;20:391–9. [4] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniquestheory and applications in engineering. Berlin, Heidelberg, New York: Springer; 1984. [5] Ioakimidis NI. Application of finite-part integrals to the singular integral equations of crack problems in plane and three-dimensional elasticity. Acta Mech 1982;45:31–47. [6] Tsamasphyros G. Numerical methods for fracture parameters calculation. In: Carpinteri A, editor. Handbook of fatigue crack propagation in metallic structures. Amsterdam: Elsevier Science B.V.; 1994. [7] Theotokoglou EE, Tsamasphyros G. Exact expressions for twodimensional singular and finite part integrals applicable in solid mechanics. Eng Anal Bound Elem 1998;22:125–32. [8] Mukherjee S. CPV and HFP integrals and their applications in the boundary element method. Int J Solids Struct 2000;37:6623–34. [9] Cruse TA, Aithal R. Non-singular boundary integral equation implementation. Int J Numer Methods Eng 1993;36:237–54. [10] Krishnasamy G, Rizzo FJ, Liu YJ. Boundary integral equations for thin bodies. Int J Numer Methods Eng 1994;37:107–21. [11] Tsamasphyros G, Theocaris PS. Integral-equation solution for half planes bonded together or in contact and containing internal cracks or holes. Ing Arch 1983;53:225–41. [12] Theotokoglou EN, Tsamasphyros G. Integral equations for any configuration of curved cracks and holes in an elastic strip. Ing Arch 1987;57:3–15. [13] Dirgantara T, Aliabadi MH. Crack growth analysis of plates loaded by bending and tension using dual boundary element method. Int J Fract 2000;105:27–74. [14] Aliabadi MH, Martin D. Boundary element hyper-singular formulation for elasto-plastic contact problems. Int J Numer Methods Eng 2000;48:995–1014. [15] Muskhelishvili NI. Singular integral equations. Groningen, The Netherlands: P. Noordhoff Ltd.; 1953. [16] Sladek N, Sladek J, Tanaka M. Regularization of hypersingular and nearly singular integrals in the potential theory and elasticity. Int J Numer Methods Eng 1993;36:1609–28. [17] Liu YJ, Rudolphi TJ. New identities for fundamental solutions and their applications to non-singular boundary element formulations. Comput Mech 1999;24:286–92.
ARTICLE IN PRESS 766
E.E. Theotokoglou, G. Tsamasphyros / Engineering Analysis with Boundary Elements 30 (2006) 758–766
[18] Mukherjee S. Finite parts of singular and hypersingular integrals with irregular boundary source points. Eng Anal Bound Elem 2000;24: 767–76. [19] Mukherjee S, Chati MK, Shi XL. Evaluation of nearly singular integrals in boundary element contour and node methods for threedimensional linear elasticity. Int J Solids Struct 2000;37:7633–54. [20] Telles JCF. A self-adaptive co-ordinate transformation for efficient numerical evaluation of general boundary element integrals. Int J Numer Methods Eng 1987;24:959–73. [21] Krishnasamy G, Schmerr LW, Rudolphi TJ, Rizzo FJ. Hypersingular boundary integral equations: some applications in acoustic and elastic waves scattering. ASME J Appl Mech 1990;57:404–14. [22] Dumont NA. On the efficient numerical evaluation of integrals with complex singularity poles. Eng Anal Bound Elem 1994;13:155–68. [23] Yun BI. A composite transformation for numerical integration of singular integrals in the BEM. Int J Numer Methods Eng 2003;57: 1883–98. [24] Kim P, Choi UJ. Two trigonometric quadrature formulas for evaluating hypersingular integrals. Int J Numer Methods Eng 2003;56:469–86. [25] Huang Q, Cruse TA. Some notes on singular integral techniques in Boundary Element Analysis. Int J Numer Methods Eng 1993;36: 2643–51. [26] Johnston PR. Application of sigmoidal transformations to weakly singular and nearly-singular boundary element integrals. Int J Numer Methods Eng 1999;45:1333–48. [27] Johnston PR, Elliott D. Error estimation of quadrature rules for evaluating singular integrals in boundary element problems. Int J Numer Methods Eng 2000;48:949–62. [28] Ma H, Kamiya N. A general algorithm for the numerical evaluation of nearly singular boundary integrals of various orders for two- and three-dimensional elasticity. Comput Mech 2002;29:277–88. [29] Korneichuk AA. Quadrature formulae for singular integrals. Z Vychisl Mat I Mat Fiz 1964;4/4(suppl.):64–74 [in Russian]. [30] Paget DE, Elliott D. An algorithm for the numerical evaluation of certain Cauchy principal value integrals. Numer Math 1972;19: 373–85. [31] Hunter DB. Some Gauss-type formulae for the evaluation of Cauchy principal values of integrals. Numer Math 1972;19:419–24. [32] Chawla MM, Ramakrishnan TR. Modified Gauss–Jacobi quadrature formulas for the numerical evaluation of Cauchy type singular integrals. BIT 1974;14:14–21. [33] Tsamasphyros G, Theocaris PS. Numerical solution of systems of singular integral equations with variable coefficients. Report.
[34]
[35]
[36] [37]
[38]
[39]
[40] [41]
[42] [43] [44] [45] [46] [47] [48] [49] [50]
National Technical University of Athens, November 1976, also. Appl Anal 1979;9:37–52. Tsamasphyros G, Theocaris PS. Sur une methode generale de quadrature pour des integrals du type Cauchy. Symp. Appl. Mathem., Salonica August 1976, also (1980). Rev Roum Sci Techn, Ser Mech Appl 1980;25:839–56. Ioakimidis NI, Theocaris PS. On the numerical evaluation of Cauchy principal value integrals. Symp. Of Appl. Mathem., Salonica August 1976, also (1977). Rev Roum Sci Techn, Ser Mech Appl 1977;22:803–18. Monegato G. On polynomials orthogonal with respect to particular variable-signed weight functions. ZAMP 1980;31:459–555. Ninis T. Factors influencing the solution of singular integral equations, Master Thesis (in Greek), National Technical University of Athens, 1980. Strouboulis T. Problems in the solution of singular integral equations, Master thesis (in Greek), National Technical University of Athens, 1980. Ioakimidis NI. A direct mehod for the construction of Gaussian quadrature rules for Cauchy type and finite part integrals. L’ Anal Numer Theorie Approx 1983;12:131–40. Monegato G. Quadrature formulas for functions with poles near the interval of integration. Math Comput 1986;47:137–58. Tsamasphyros G. A study of factors influencing the solution of singular integral equations. Eng Fract Mech 1986;24: 567–78. Tsamasphyros G, Dimou G. A Gauss quadrature formula for Cauchy type integrals. Comput Mech 1989;4:137–48. Tsamasphyros G, Dimou G. Gauss quadrature rules for finite part integrals. Int J Numer Methods Eng 1990;30:13–26. Laurie DP. Computation of Gauss-type quadrature formulas. J Comput Appl Math 2001;127:201–17. Isaacson E, Keller HB. Analysis of numerical methods. London: Wiley; 1966. Hildebrand FB. Introduction to numerical analysis. 2nd ed. New York: McGraw-Hill; 1974. Szego G. Orthogonal polynomials. New York: American Mathematical Society; 1959. Guelfond AO. Calcul des differences finies. Paris: Dunod; 1963. Stroud A, Secret D. Gaussian quadrature formulas. Englewood Cliffs, N.J: Prentice-Hall; 1966. Gladwell GML, England AH. Orthogonal polynomial solutions to some mixed boundary value problems in elasticity theory. Q J Mech Appl Math 1977;30:175–85.