• 検索結果がありません。

# Numerical Integration of Functions with Endpoint Singularities and/or Complex Poles

N/A
N/A
Protected

シェア "Numerical Integration of Functions with Endpoint Singularities and/or Complex Poles"

Copied!
27
0
0

(1)

### in 3D Galerkin Boundary Element Methods

To Professor Masatake Mori for his outstanding contribution to numerical integration

By

Giovanni Monegato and Letizia Scuderi∗∗

Abstract

In this paper we propose special strategies to compute 1D integrals of functions having weakly or strong singularities at the endpoints of the interval of integration or complex poles close to the domain of integration. As application of the proposed strategies, we compute a four dimensional integral arising from 3D Galerkin boundary element methods (BEM) applied to hypersingular boundary integral equations.

§1. Introduction

In the computation of integrals and in the numerical solution of integral equations, one often has to deal with the numerical integration of functions with endpoint weak singularities, or with analytic functions having complex conjugate poles near the domain of integration. Nonlinear changes of variables have shown to be a very eﬀective tool to face these problems. Indeed, in the case of weak singularities a proper change of variable can make the integrand function arbitrarily smooth. This, combined with the use of a standard rule

Communicated by H. Okamoto. Received September 28, 2004. Revised February 22, 2005.

2000 Mathematics Subject Classiﬁcation(s): 65D32, 65R20, 65N38.

This work was supported by the Ministero dell’Istruzione, dell’Universit`a e della Ricerca of Italy and by the INDAM 2003 project “Modellistica Numerica per il Calcolo Scientiﬁco e Applicazioni Avanzate”.

∗∗Dipartimento di Matematica, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy.

(2)

like the Gauss-Legendre one or the composite trapezoidal formula, allows to achieve high accuracy using a low number of abscissas (see [2], [8], [9], [11], [14], [16]). Here we recall some of the known changes of variable, all having as co-domain the interval (0,1), which have been proposed in the literature. For some comparisons (see [9]).

In [8], the following polynomial transformation, which is a generalization of the one proposed in [6] to construct lattice rules for multiple integration on the hypercube, has been used:

(1.1) ϕ1(t) = (p+q−1) (p1)(q1)

t 0

up1 (1−u)q1du, 0≤t≤1, p, q1.

The integral in (1.1) could be computed, for any t (0,1], using recurrence relations (see [9]). In any case it can be evaluated more eﬃciently, for the given values of p and q, by means of the n-point Gauss-Legendre rule, with n=p+q2 .

Then, we recall a trigonometric transformation of the form (1.2) ϕ2(t) =

t

0(sinπ2u)p1(cosπ2u)q1du 1

0(sinπ2u)p1(cosπ2u)q1du, 0≤t≤1, p, q 1.

This latter generalizes a trigonometric transformation proposed in [14]. Also in this case the integral in (1.2) could be computed for any t∈(0,1] by using proper recurrence relations or by a direct computation using a Gauss-Legendre rule.

A third transformation, of rational type, used in [11], [3] in connection with the solution of certain integral equations, and in [4] to evaluate Hadamard ﬁnite-part integrals, is the following one:

(1.3) ϕ3(t) = tp

tp+ (1−t)q, 0≤t≤1, p, q1.

The last transformation we recall is of exponential type:

(1.4) ϕ4(t) = 1 2+1

2tanh π

2sinh(t)

, −∞< t <∞.

It is the well-known double exponential (DE) transformation (see [16], [10], [15]).

In this paper we examine the use of nonlinear changes of variable for the numerical evaluation of some 4-dimensional integrals arising in the numerical solution of hypersingular boundary integral equations. In particular we consider

(3)

3D Galerkin boundary element methods, where a key issue is the eﬃcient and fast computation of integrals of the form

(1.5)

i

j

k(x, y)

||x−y||3dy dx,

where ∆i, ∆jare triangular elements, possibly curved, andk(x, y) is bounded.

When ∆i and ∆j are disjoint and not too close to each other, the compu- tation of the above integral can be performed by standard rules. The diﬃcult cases are those where (i) ∆i and ∆j are coincident, or share a common edge, or (ii) are disjoint but very close to each other. When ∆i and ∆j have only a common point, but are not close to each other at some other points, the computation can be classiﬁed standard.

Many papers have been devoted to this topic (see, for example, [1], [5], [12], [13]), although none of them seems to treat the case (ii), which is considered by some users the most “diﬃcult” case. However in our opinion none of them are fully satisfactory because generally they either require some analytical calcula- tion which is by non means trivial, or do not take into account the presence of complex poles which may eﬀect adversely the numerical computation.

Indeed the strategy adopted by these methods is that of reducing the computation of the multidimensional integral to the evaluation of nested 1D integrals of analytic functions. However this is still not satisfactory because it is by now well-known that when integrating an analytic function, the position of its poles is crucial for the performance of a quadrature rule. When an analytic function has poles which are outside the domain of integration but lie very close to it, and this is what happens in the computation of some of our integrals, its numerical integration by means of a quadrature rule can give very poor accuracy if one does not adopt a special strategy.

In this paper we examine the case of (1.5) with ∆i, ∆j coinciding. The other cases of ∆i, ∆j having a common edge or disjoint but very close to each other will be treated in a following paper. We believe that our strategy can be applied successfully also to these cases.

§2. One-dimensional Integration Rules

There are two typical approaches to compute (1.5). One uses the global cartesian coordinates (see [5]), the other uses local polar coordinates (see [12]).

Both of them reduce (1.5) to four nested one-dimensional integrals, whose in- tegrand functions are analytical in the ﬁrst case, and analytical or with log- singularities at the endpoints of the interval of integration in the second case.

(4)

Unfortunately, in both cases the generated analytical functions have complex poles, depending upon outer variables or the form of the triangles, which may be arbitrarily close to the interval of integration. Therefore, for both approaches the Gaussian quadrature rules perform very poorly unless one considers a high number of quadrature nodes. Recalling that we have to compute four nested one-dimensional integrals, we can use only a few (let us say, at the most 6÷10) quadrature nodes for each integral. Notice that 10 quadrature nodes for each integral means a total of 104quadrature nodes to compute (1.5), which repre- sents only one element of the system matrix arising from the Galerkin boundary element method. For this reason in this section we propose special strategies, which allow us to compute the integrals with the above-mentioned character- istics using very few Gaussian quadrature nodes. In particular we propose to use variable transformations of the type (1.1), either to make the integrand functions smoother at the endpoints of the interval of integration, or to move poles away from the interval. Then, the Gauss-Legendre rule is applied.

In these situations the DE-rule, which is obtained by combining (1.4) with the (truncated) trapezoidal rule (see [16]), has certainly a much higher (ex- ponential) rate of convergence. However, in the particular application we will consider in Section 3, it is of key importance to obtain a good accuracy with a very small number of function evaluations. For such small values of n, let us sayn≤10, our approach appears more eﬃcient.

A typical case is an integral of a function having only a log-singularities at the endpoints of the interval of integration. For example,

(2.1) I1=

1 0

exlog1−x x dx.

By introducing in (2.1) the change of variable (1.1), that here we write in the form

(2.2) x=γ0,1q0,q1(t) := (q0+q11) (q01)(q11)

t 0

xq01(1−x)q11dx, q0, q11, and then applying then-point Gauss-Legendre rule we obtain the improvements reported in Table 1. Notice that in (2.2) forq1= 1 we have

(2.3) γ0,1q0,1(t) =tq0 =:γ0q0(t) and forq0= 1

(2.4) γ0,11,q1(t) = 1(1−t)q1 =:γ1q1(t);

(5)

the choices (2.3) and (2.4) allow to smooth the singularities at the endpoints 0 and 1, respectively. Ifq0=q1=q, we set

(2.5) γ0,1q :=γ0,1q,q.

n q= 1 q= 2 q= 3 q= 4

2 2.5001 5.6501 9.0701 9.0401 4 6.6702 1.3502 9.5603 1.3401 8 1.8002 1.2103 1.8704 6.8305 16 4.7203 8.3205 3.0906 2.0007 32 1.2103 5.4906 5.1708 8.3410 64 3.0804 3.5307 8.4010 3.3612 128 7.7505 2.2408 1.3111 2.2713 256 1.9505 1.4109 7.8114 1.3813

512 4.8706 8.8911 −− −−

Table 1. Relative errors inI1=1

0 eγ(t)log1γ(t)γ(t)γ(t)dt,γ(t) =γ0,1q (t).

In all tables the sign “−−” means that full relative accuracy (i.e. 14 signiﬁcant digits in our case) has been achieved.

Another situation that we will meet is that of integrals of the form

(2.6) I2=

1 0

f(x) x2+ε2dx

withεsmall andf(x) smooth. To compute eﬃciently this integral it is conve- nient to introduce preliminarily the simple change of variable

x=γ0q(t), q >1,

and then to apply the n-point Gauss-Legendre rule. For example, if we take f(x) =ex we obtain the relative errors of Table 2.

ε= 101 ε= 103 ε= 105

n q= 1 q= 2 q= 1 q= 4 q= 1 q= 7

4 4.2502 6.1103 9.7201 6.3101 1.00 + 00 8.8801 8 1.4503 5.7504 9.0501 1.7101 9.9901 2.7401 16 9.7507 1.2007 6.6101 2.7303 9.9601 1.1101 32 1.9813 −− 7.2002 3.1605 9.8701 2.8503 64 −− −− 6.3003 1.2910 9.4701 4.5007

128 −− −− 3.5305 −− 7.9201 1.4312

256 −− −− 3.4410 −− 2.9001 −−

512 −− −− −− −− 3.8502 −−

Table 2. Relative errors inI2 =1 0

qtq−1etq t2q2 dt.

Since our 1D rules will be applied to a quadruple integral, it is important to obtain the required accuracy using a small number of nodes. In the integral

(6)

(2.6), whenε is very small, for example in the above caseε= 105, the value ofqone has to take makes the new integrand function so ﬂat around the origin, that to achieve a good accuracy one has to take a value ofnwhich is not any longer small. In this situation it is more convenient to proceed as follows. Write

(2.7) I2=

ε 0

+ 1

ε

f(x) x2+ε2dx;

then apply the n-point Gauss-Legendre rule to the ﬁrst integral on the right- hand side, and to the second one after having introduced in it the change of variable x = γ0q(t). Subdivision (2.7) is of key importance for the following reasons. In (2.6) the transformation x=γ0q(t) leaves unchanged the interval of integration, but as q gets larger, the modulus of the poles tends to 1, and after a certain value of q some of the poles move again towards the interval of integration. In (2.7) this unwanted phenomenon does not arise. Indeed, in the ﬁrst integral on the right-hand side the distance of the two poles from the interval of integration is equal to the length of this latter, while the length of the interval of integration of the second integral, after having introduced the transformationx=γ0q(t), tends to zero asqtends to inﬁnity. The eﬀectiveness of this new approach is shown by the results in Table 3.

ε= 10−1 ε= 10−3 ε= 10−5

2n q= 3 q= 100 q= 4 q= 100 q= 7 q= 100

4 8.2603 6.0403 3.0001 1.2802 4.4001 2.0401 8 6.0705 1.1505 1.3103 4.5003 8.3602 3.4202 16 5.4409 2.1109 2.9304 3.2005 9.4603 4.9004 32 −− −− 7.7407 5.0110 2.3805 7.6908

64 −− −− 6.5413 −− 5.0510 −−

Table 3. Relative errors inI2=ε 0

ex

x22dx+1 ε 1 q

qtq−1etq t2q2 dt.

Similarly, if the integral has the form

(2.8) I3=

b a

f(x) (x2+ε2)αdx

witha <0,b >0 and α >0, then we proceed as in (2.7) by splitting (2.8) as follows

(2.9) I3=

ε a

+ ε

ε

+ b

ε

f(x) (x2+ε2)αdx.

We apply the n-point Gauss-Legendre quadrature rule to each integral, after having introduced the change of variable x= tq in the ﬁrst and in the third

(7)

integral to the right-hand side of (2.9). In the following Table 4, we have reported the relative errors obtained by this last numerical procedure when in (2.8)f 1,a=0.5,b= 0.5 andα= 3/2.

In Table 4 and belowGLmdenotes them-point Gauss-Legendre rule.

ε= 10−1 ε= 10−3 ε= 10−5

3n GL3n q= 2 GL3n q= 4 q= 100 GL3n q= 7 q= 100 6 3.1501 5.9802 1.00 + 00 3.1101 1.5601 1.00 + 00 3.4801 3.0401 12 4.1702 2.3103 9.9901 2.9602 8.0303 1.00 + 00 1.5701 1.3602 18 4.6203 8.2705 9.9901 9.1603 7.7204 1.00 + 00 1.8703 4.9803 24 4.8504 2.7706 9.9801 6.7504 3.0205 1.00 + 00 1.3002 9.4404 30 4.9605 9.0608 9.9701 2.2004 9.7007 1.00 + 00 1.7403 1.9505 36 4.9706 2.9109 9.9501 1.3405 3.0208 1.00 + 00 6.0804 1.3705 72 4.2412 −− 9.8201 6.1211 −− 1.00 + 00 1.3807 5.8211

144 −− −− 9.3201 −− −− 1.00 + 00 −− −−

Table 4. Relative errors inI3= (|a|1q

ε 1 q

+b 1q ε

1

q)qtq1(t2q+ε2)3/2dt +ε

ε(x2+ε2)−3/2dx,a=0.5,b= 0.5.

The integral we have to compute could also be of the type

(2.10) I4=

1 0

f(x)

[(x−r)2+ε2]αdx

with 0 < r <1 andα > 0. Notice that also I3 can be rewritten in this form and vice versa. In this case we can proceed as follows. Write ﬁrst

I4= r

0

+ 1

r

f(x)

[(x−r)2+ε2]αdx, and then

(2.11) I4= 1−r r

1 r

f(1rr(1−t))

[(t−r)2+ (1rrε)2]αdt+ 1

r

f(x)

[(x−r)2+ε2]αdx.

This splitting is suggested by the remark that, when one uses Gaussian rules, poles with real parts close to the endpoints of the interval of integration are less adverse than poles with a real part in the central part of the same interval.

Finally, we introduce, as forI2, the change of variablet, x=uq,q≥1, in both integrals and apply the same Gauss-Legendre rule to each of them.

Since this strategy is also alternative to a splitting of the type (2.9), we have applied it to the same integral considered in Table 4. Therefore, in Table 5 we have reported the relative errors obtained by taking in (2.10)f 1,r= 0.5 andα= 3/2.

(8)

ε= 10−1 ε= 5·10−2

2n GL2n q= 1 q= 50 GL2n q= 1 q= 50

6 3.1501 6.7701 4.9502 7.1401 1.4601 3.8002 12 4.1702 2.4903 2.5506 3.3401 2.2502 3.0603 18 4.6203 6.0505 9.7006 1.2701 2.1903 1.8304 24 4.8504 1.0306 4.4808 4.4302 1.5404 1.0005 30 4.9605 9.0809 1.2809 1.4802 7.7206 5.1807 36 4.9706 1.6010 1.2111 4.8603 1.6207 2.5808

72 4.2412 −− −− 5.0706 −− −−

144 −− −− −− 4.0412 −− −−

Table 5. Relative errors inI4=1rr1 r

1

qquq1f(1−rr (1−uq))

(uq−r)2 +(1−rr ε)23/2

du+1 r

1

qquq−1f(uq)[(uq−r)2+ε2]−3/2du, r= 0.5.

From a comparison between the Tables 4 and 5, it appears that (2.11) is more eﬀective when ε is not too small; let us say of order greater than 102. Otherwise (2.9) seems to produce better results.

The above strategy has been applied also to similar integrals with the algebraic polynomial (of degree 2) replaced by a corresponding trigonometric polynomial. In particular, we have also considered integrals of the type

(2.12) I5= b

a

f(x)

[(a1cos(x) +b1sin(x))2+ (a2cos(x) +b2sin(x))2]αdx.

Byz0we denote the pair of complex conjugate poles of the trigonometric poly- nomial closest to the interval of integration. In this situation, we have applied the above strategy by choosing r=Re(z0). The results we have obtained by taking in (2.12) f 1, a= 1.5, b = 0, α= 3/2 and a1 = 0, a2 =b2 = 1, b1=c with varyingc, are reported in Table 6.

c= 1 c= 1/2 c= 1/4 c= 1/8

Re(z0)≈ −0.55 Re(z0)≈ −0.72 Re(z0)≈ −0.77 Re(z0)≈ −0.78 Im(z0)≈ ±0.40 Im(z0)≈ ±0.24 Im(z0)≈ ±0.12 Im(z0)≈ ±0.06 2n GL2n q= 50 GL2n q= 50 GL2n q= 50 GL2n q= 50

8 6.4704 3.1904 2.4202 2.4903 2.4301 8.9303 6.0301 8.4602 16 4.7008 1.1607 1.1204 4.8307 1.7802 1.1004 7.8602 2.7703 32 1.5714 1.9114 9.5509 7.5413 2.3705 9.8210 2.5602 1.6307

64 −− −− −− −− 7.1509 −− 1.1904 3.1714

128 −− −− −− −− −− −− 1.4609 −−

Table 6. Relative errors inI5=0

1.5[(csin(x))2+ (cos(x) + sin(x))2]3/2dx.

(9)

Also integrals of the form

(2.13) I6=

1 ε

f(x) x dx,

withf(x) smooth andε >0 small, can be eﬃciently evaluated using a Gauss- Legendre rule, after having introduced the change of variable x = tq. Some relative errors concerning an integral of this type, withf(x) =ex, are reported in Table 7.

ε= 10−1 ε= 10−3 ε= 10−5

n q= 1 q= 3 q= 1 q= 6 q= 1 q= 8

4 4.4402 2.3606 3.3801 2.3704 5.7301 1.0303 8 2.4805 2.7512 1.9601 8.4709 4.7401 2.9407

16 7.1510 −− 7.4502 −− 3.7101 −−

32 −− −− 1.0902 −− 2.6601 −−

64 −− −− 2.0504 −− 1.6401 −−

128 −− −− 6.4308 −− 7.3702 −−

256 −− −− −− −− 1.6002 −−

Table 7. Relative errors inI6=1

ε 1 q

qetq t dt.

In the introduction we have remarked that for small values ofn, the ap- proach we have proposed gives slightly better results than the DE-rule. This is due to the fact that the stepsize hof the DE-rule is chosen to optimize the rate of convergence, that is, the behavior of the remainder term asn→ ∞(see [15]). Indeed, asn→ ∞the accuracy given by the DE-rule is superior to that given by our numerical integration approach.

However, in the applications we are interested in, we ought to achieve a good accuracy by taking a small value ofn. Therefore, in the case of the DE- rule, for a given (low) value ofnwe would like to choose the parameterhwhich minimizes the error.

For the integralsI1throughI5of this section we have plotted the behavior of the error term as a function ofh, for several values ofn. In all cases we have examined, a proper choice of hhas reduced substantially the error. Here we report a few examples, which are however representative of all cases we have examined. In all cases there exists an optimal choice ofhwhich allows to reach very high accuracies using a small value ofn. However we do not know a simple formula to compute such value, or a good approximation of it. We note that all ﬁgures have been obtained by subdividing the domain ofh into 100 equal parts. Therefore the actual minima could be even smaller than those appearing in the ﬁgures.

(10)

0.45 0.5 0.55 0.6 0.65 10−4

10−3 10−2

stepsize h

|relative error|

Figure 1. IntegralI1,n= 3,h∈(0.45,0.65).

0.35 0.4 0.45 0.5 0.55

10−7 10−6 10−5 10−4

stepsize h

|relative error|

Figure 2. IntegralI1,n= 5,h∈(0.35,0.55).

(11)

4.72 4.74 4.76 4.78 4.8 4.82 4.84 4.86 4.88 4.9 x 10−3 10−10

10−9 10−8 10−7

stepsize k (h=k+0.48)

|relative error|

Figure 3. IntegralI1,n= 5,h∈(0.4847,0.4849).

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

10−4 10−3 10−2 10−1

stepsize h

|relative error|

Figure 4. IntegralI4,r= 0.5,ε= 101,n= 3,h∈(0.1,2).

(12)

2.89 2.891 2.892 2.893 2.894 2.895 2.896 2.897 2.898 2.899 2.9 x 10−5 10−10

10−9 10−8

stepsize k (h=k+1.1355)

|relative error|

Figure 5. IntegralI4,r= 0.5,ε= 101,n= 3,h∈(1.1355289,1.135529).

2.38 2.382 2.384 2.386 2.388 2.39 2.392 2.394 2.396 2.398 2.4

x 10−6 10−9

10−8 10−7 10−6

stepsize k (h=k+1.48e−03)

|relative error|

Figure 6. IntegralI4,r= 0.5,ε= 103,n= 3,h∈(1.48238e03,1.4824e03).

(13)

6.35 6.355 6.36 6.365 6.37 6.375 6.38 6.385 6.39 6.395 6.4 x 10−5 10−11

10−10 10−9 10−8

stepsize k (h=k+3.9797)

|relative error|

Figure 7. IntegralI5,ε= 1/2,n= 3,h∈(3.9797635,3.979764).

8.32 8.321 8.322 8.323 8.324 8.325 8.326 8.327 8.328 8.329 x 10−3 10−9

10−8 10−7 10−6

stepsize k (h=k+2.40)

|relative error|

Figure 8. IntegralI5,ε= 1/8,n= 3,h∈(2.40832,2.40833).

(14)

The last integral we need to consider is

(2.14) I7=

R(y) 0

= f(x) x dx,

deﬁned as Hadamard ﬁnite-part, wherey is an outer variable, which can make R(y) very close to zero. Since a change of variable does not leave unchanged the value of the integral, and we do not want to construct an ad hoc rule for each value of the variable y, we suggest to proceed as follows. If f(0) exists and can be easily computed, then write

(2.15)

R(y) 0

= f(x) x dx=

R(y) 0

f(x)−f(0)

x dx+f(0) logR(y)

and apply the Gauss-Legendre to the integral on the right-hand side of (2.15).

If however this is not the case, then write

(2.16)

R(y) 0

= f(x) x dx=

r 0

= +

R(y) r

f(x)

x dx.

where r > 0 has a suitable ﬁxed value. To compute the ﬁrst integral on the right-hand side we use a standard product rule for Hadamard ﬁnite part in- tegrals, whose coeﬃcients have been computed once for all ([7]). The sec- ond integral can be computed by the technique suggested for (2.13) with a large q (for example, q = 100) when f is smooth. The choice of r may depend on f. Indeed, if f has a pair of complex poles with the real part ﬁxed, independent ofy and belonging to the interval of integration (0, R(y)), it is convenient to take r equal to the real part of the poles. If the real part varies with y in a little range of (0, R(y)) we can choose r as an in- termediate value of the range; otherwise, it is convenient to use a product formula without splitting. Finally, if f is smooth, we suggest to choose r small.

In the Tables 8 and 9 we have reported the relative errors obtained ﬁrst by applying directly the product rule with n quadrature points (P Rn) to I5 in (2.14), and then by splitting I5 as in (2.16) and computing the in- tegrals to the left-hand side as suggested above (P Rn +GLn). In partic- ular, in Table 8 we report the relative errors corresponding to the choice f(x) = f1(x) = 1/[(x−r)2 +ε2] with r = 0.1 and ε varying; in Table 9 we give those corresponding tof(x) =f2(x) =ex andr= 0.1 in (2.16).

(15)

ε= 10−1 ε= 10−2 ε= 10−3 2n P R2n P Rn+GLn P R2n P Rn+GLn P R2n P Rn+GLn

4 1.83 + 00 2.3101 1.00 + 00 4.7201 1.00 + 00 8.3501 8 2.45 + 00 1.5001 3.15 + 00 1.03 + 00 1.27 + 00 4.1701 16 9.6202 2.8604 1.06 + 01 6.5402 6.66 + 00 9.0101 32 6.5305 9.9210 6.86 + 00 7.6003 9.63 + 00 1.07 + 00 64 1.3011 2.3213 1.66 + 00 1.5105 9.03 + 00 1.6401

128 −− −− 1.4201 3.3312 3.67 + 01 2.3103

256 −− −− 4.5404 −− 4.45 + 01 5.7707

512 −− −− 9.0310 −− 1.22 + 00 3.8615

Table 8. Relative errors inI7=R(y) 0

[(x0.1)22]−1

x dx

= (r

0 +R(y)

r )[(x0.1)x22]−1dxwithr= 0.1,R(y) = 0.5.

R(y) = 10−2 R(y) = 0.5

2n P R2n P Rn+GLn P R2n P Rn+GLn

4 1.1411 1.1303 1.8003 1.9702 8 −− 1.2907 4.6910 6.4906 16 −− 6.3815 4.0913 1.2012

32 −− −− −− 3.6513

64 −− −− −− −−

Table 9. Relative errors inI7=R(y)

0 ex

x dx=r

0 ex

x dx +(R(y))

1 q r

1q

qetq

t dtwithr= 0.1,q= 100.

§3. The Four-dimensional Integral

The integral we have to compute has the following form

(3.1) I=

= d¯x

=Nix)Njy)

||x¯−y¯||3 d¯y,

where for simplicity we assume ∆ to be a triangle with corners at (0,0), (a1, a2), (b1, b2) and{Nix)},{Njy)}are the shape and the test functions, respectively.

Our approach however easily extends to the case of a curved triangle, given by a smooth parametric representation.

By the linear transformation ¯x=Ax, ¯y=Aywith A=

a1 b1

a2 b2

we map the domain ∆ into the reference triangle T having corners at (0,0), (1,0) (0,1). In this case the integral (3.1) can be replaced by

(16)

(3.2) IT =|A|2

T

= dx

T

= Ni(x)Nj(y)

||A(x−y)||3dy, where|A|= det(A),Ni(x) :=Ni(Ax) andNj(y) :=Nj(Ay).

The value ofIT is generally diﬀerent from that ofI; however the sum of the values of the integrals deﬁned over all boundary elements ∆i is equal to the sum of the corresponding integrals of form (3.2). Therefore from now on we refer to the form (3.2).

To compute (3.2) we split it as follows:

IT =|A|2

T

Ni(x)dx

T

−Nj(y)−Nj(x)

||A(x−y)||3 dy (3.3)

+

T

Ni(x)Nj(x)dx

T

= dy

||A(x−y)||3

=:|A|2[IT ,1+IT ,2]

since the inner integral in IT ,2 can be evaluated analytically and the strong singularity inIT ,1is of Cauchy type.

As regardsIT ,1, ﬁrst we introduce the following polar coordinates (3.4)

y1=x1+rcos(ϑ) y2=x2+rsin(ϑ)

where 0≤r≤R(ϑ), 0≤ϑ≤2π. Then we split the interval (0,2π) into three subintervals (ϑ0, ϑ1), (ϑ1, ϑ2) and (ϑ2, ϑ3), with

(3.5)



























ϑ0=−π+ arctan x2

x1

ϑ1=arctan x2

1−x1

ϑ2=π−arctan

1−x2 x1

ϑ3=π+ arctan x2

x1

.

Forϑ∈l1, ϑl],l= 1,2,3, we have 0≤r≤Rl with (3.6) Rl:=Rl(x1, x2, ϑ) =dl(x1, x2)

pl(ϑ)

(17)

where (3.7)









d1(x1, x2) =x2

d2(x1, x2) = 1

2(1−x1−x2) d3(x1, x2) =x1

,









p1(ϑ) =sin(ϑ) p2(ϑ) = 1

2[sin(ϑ) + cos(ϑ)]

p3(ϑ) =cos(ϑ)

.

Our integralIT ,1 can then be expressed as follows:

IT ,1= 1

0

dx1

1x1 0

Ni(x1, x2)dx2

(3.8)

3 l=1

ϑl

ϑl−1

Rl

0

= Nj(x1+rcos(ϑ), x2+rsin(ϑ))−Nj(x1, x2)

r2f(ϑ) dr,

where

f(ϑ) = [(a1cos(ϑ) +b1sin(ϑ))2+ (a2cos(ϑ) +b2sin(ϑ))2]3/2. Notice thatf(ϑ) vanishes when

tan(ϑ) =α±iβ, with

α=−a1b1+a2b2

b21+b22 , β =a2b1−a1b2

b21+b22 .

In particular the poles off(ϑ) depend on the corners of the triangle ∆. There- fore their distance from the interval of integration with respect to the variable ϑ depends upon the shape of ∆ and could be very small. For example, in the case of a triangle ∆, whose corners are at (0,0), (a,1), (0,1), we have ap- proximately the following complex poles: (1.01 +kπ)±0.40ı when a = 1, (0.85 +kπ)±0.24ı whena= 1/2 and (0.80 +kπ)±0.12ı whena= 1/4, k∈Z .Z

To proceed we set

Njr(x1, x2;r, ϑ) = Nj(x1+rcos(ϑ), x2+rsin(ϑ))−Nj(x1, x2) r

(18)

and write (3.9)

IT ,1= 1

0

dx1 1x1

0

Ni(x1, x2)dx2 3

l=1

ϑl ϑl−1

[f(ϑ)]1 Rl

0

= Njr(x1, x2;r, ϑ)

r dr

= 1

0

dx1 1x1

0

Ni(x1, x2)dx2 3 l=1

ϑl

ϑl−1

[f(ϑ)]1 Rl

0

Njr(x1, x2;r, ϑ)−Njr(x1, x2; 0, ϑ)

r dr+Njr(x1, x2; 0, ϑ) log(Rl)

Remark. Notice that in the case of {Nj} linear, that is Nj =Njl, j = 1,2,3, and

(3.10)

N1l(x1, x2) =x1, N2l(x1, x2) =x2,

N3l(x1, x2) = 1−x1−x2, we have

N1r(x1, x2;r, ϑ) = cos(ϑ), N2r(x1, x2;r, ϑ) = sin(ϑ),

N3r(x1, x2;r, ϑ) =−[sin(ϑ) + cos(ϑ)].

Therefore Njr depends only upon ϑ and the integral over (0, Rl) in the ﬁrst expression ofIT ,1 in (3.9) can be performed analytically.

By denoting

Dj(x1, x2;r, ϑ) :=Njr(x1, x2;r, ϑ)−Njr(x1, x2; 0, ϑ)

r ,

we write

IT ,1= 3 l=1

IT ,1l ,

The first case is the Whitham equation, where numerical evidence points to the conclusion that the main bifurcation branch features three distinct points of interest, namely a

This review is devoted to the optimal with respect to accuracy algorithms of the calculation of singular integrals with ﬁxed singu- larity, Cauchy and Hilbert kernels, polysingular

This review is devoted to the optimal with respect to accuracy algorithms of the calculation of singular integrals with ﬁxed singu- larity, Cauchy and Hilbert kernels, polysingular

This review is devoted to the optimal with respect to accuracy algorithms of the calculation of singular integrals with ﬁxed singu- larity, Cauchy and Hilbert kernels, polysingular

Sofonea, Variational and numerical analysis of a quasistatic viscoelastic problem with normal compliance, friction and damage,

Then it follows immediately from a suitable version of “Hensel’s Lemma” [cf., e.g., the argument of [4], Lemma 2.1] that S may be obtained, as the notation suggests, as the m A

A mathematical formulation of well-posed initial boundary value problems for viscous incompressible fluid flow-through-bounded domain is described for the case where the values

We provide an efficient formula for the colored Jones function of the simplest hyperbolic non-2-bridge knot, and using this formula, we provide numerical evidence for the