Numerical integrations over an arbitrary ... - Semantic Scholar
Applied Mathematics and Computation 210 (2009) 515–524
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Numerical integrations over an arbitrary quadrilateral region Md. Shafiqul Islam *, M. Alamgir Hossain Department of Mathematics, University of Dhaka, Dhaka 1000, Bangladesh
a r t i c l e
i n f o
Keywords: Double integral Numerical integration Quadrilateral and triangular finite element Gaussian quadrature
a b s t r a c t In this paper, double integrals over an arbitrary quadrilateral are evaluated exploiting finite element method. The physical region is transformed into a standard quadrilateral finite element using the basis functions in local space. Then the standard quadrilateral is subdivided into two triangles, and each triangle is further discretized into 4 n2 right isosceles triangles, with area 2n12 , and thus composite numerical integration is employed. In addition, the affine transformation over each discretized triangle and the use of linearity property of integrals are applied. Finally, each isosceles triangle is transformed into a 2-square finite element to compute new n2 extended symmetric Gauss points and corresponding weight coefficients, where n is the lower order conventional Gauss Legendre quadratures. These new Gauss points and weights are used to compute the double integral. Examples are considered over an arbitrary domain, and rational and irrational integrals which can not be evaluated analytically. Ó 2009 Elsevier Inc. All rights reserved.
1. Introduction Numerical simulation in engineering science and in applied mathematics has become a powerful tool to model the physical phenomena, particularly when analytical solutions are not available and/or are very difficult to obtain. The integrals arising in practical problems are not always simple or polynomial but rational and irrational expressions, in which the quadrature scheme cannot evaluate exactly [12]. Even there is no order of Gauss quadrature that will evaluate these integrals exactly [9,10]. The integration points have to be increased in order to improve the integration accuracy and it is desirable to make these evaluations by using as few Gauss points as possible, from the point of view of the computational efficiency. Among various numerical techniques, the finite element method (FEM) is probably one of the most widely accepted even for the complex geometries. This advantage is supported by the element wise coordinate transformation from one space to the other space. From the literature review we may realize that a lot of works of numerical integration using Gauss quadrature over triangular region has been done [1–8,13], but a limited work is attempted over the quadrilateral region, such as [12]. Rathod [11] presented some analytical formulas for rational integrals over quadrilateral element but it was confined with monomials as numerators. For this, a little work has been done in this study to carry out the development of a good numerical integration technique over the arbitrary convex quadrilateral region with the advent of FEM. Very recently, a rigorous and elaborate survey has been reported in the literature [13] by Rathod et al., and they have derived various orders of extended numerical quadrature rules based on classical Gauss Legendre formula. In their work, a transformation has been used from standard triangular surface to a standard 2-square. All the formulations are derived, and all the examples are tested for certain triangular region. In contrast to this study, we use an arbitrary quadrilateral region with convex and straight sides, which is transformed into a standard square finite element by coordinate transformations. * Corresponding author. E-mail address: mdshafi[email protected] (Md. Shafiqul Islam). 0096-3003/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2009.01.030
516
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
Then the standard square is discretized into two triangular regions and each of these standard triangle is discretized into 4 n2 triangles instead of 3 n2 [13]. We then map further each of the standard triangle into the 2-square using standard quadrilateral basis functions, and taking the sum of the discretized triangles. The subsequent formulations are developed by Mathematica instead of C-Language, which can be coded easily. Examples are considered over an arbitrary domain, and a problem of particular domain which is available in the literature. The results, obtained by the present formulations, converge to the exact solutions correct upto 15 decimal places. 2. Formulation of integrals over an arbitrary quadrilateral region The integral of an arbitrary function, f(x, y) over an arbitrary quadrilateral region AQ is given by
I¼
ZZ
f ðx; yÞ dy dx ¼
ZZ
AQ
f ðx; yÞ dx dy:
ð1Þ
AQ
The integral I of Eq. (1) is then transformed into an integral over the region of the standard quadrilateral, SQ = {(u, v):1 6 u 6 1, 1 6 v 6 1}, shown in Fig. 1, by the changing the coordinates as:
x¼
4 X
xi Q i
and y ¼
i¼1
4 X
yi Q i ;
ð2aÞ
i¼1
where Qi are the bilinear quadrilateral element basis functions in (u, v)-space:
Q 1 ðu; v Þ ¼ ð1 uÞð1 v Þ=4;
Q 2 ðu; v Þ ¼ ð1 þ uÞð1 v Þ=4;
Q 3 ðu; v Þ ¼ ð1 þ uÞð1 þ v Þ=4;
Q 4 ðu; v Þ ¼ ð1 uÞð1 þ v Þ=4:
ð2bÞ
The corresponding Jacobian [11] is then
oðx; yÞ ox oy ox oy ¼ a0 þ a1 u þ a2 v ; ¼ oðu; v Þ ou ov ov ou
ð3aÞ
1 8 1 a1 ¼ ½ðx4 x3 Þðy2 y1 Þ þ ðx1 x2 Þðy4 y3 Þ; 8 1 a2 ¼ ½ðx4 x1 Þðy2 y3 Þ þ ðx3 x2 Þðy4 y1 Þ: 8
ð3bÞ
J1 ¼ where
a0 ¼ ½ðx4 x2 Þðy1 y3 Þ þ ðx3 x1 Þðy4 y2 Þ;
y 3 ( x3 , y3 )
(− 1 , 1 ) 4
v
3 ( 1, 1 )
2 (x2 , y2 )
u
( x4 , y 4 ) 4 1 ( x1 , y1 )
x
( − 1, − 1) 1
(b) Standard square element, SQ
(a) Arbitrary quadrilateral, AQ
v
(− 1 , 1 ) 4
2 (1,−1)
3 ( 1, 1 )
T '' u
T' ( − 1, − 1) 1
2 (1,−1)
(c) Discretize SQ into two triangles T ' and T '' Fig. 1. Transformation of arbitrary quadrilateral (AQ) into standard square (SQ), and discretization SQ into two triangles T0 and T00 .
517
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
v
( − 1 ,1 ) 3
1 ( 1, 1 )
v
( − 1 ,1 ) 3
T '' u
u T
2 (1, − 1)
'
2 (1, − 1)
( − 1, − 1) 1
Fig. 2. Transformation of T00 into T0 .
Now use Eqs. (1)–(3) to obtain,
I¼
ZZ
f ðx; yÞ dy dx ¼
ZZ
AQ
SQ
f ðxðu; v Þ; yðu; v ÞÞjJ 1 j du dv ¼
ZZ
gðu; v Þ du dv ¼ SQ
Z
1 1
Z
1
gðu; v Þ du dv ;
ð4Þ
1
where, g(u, v) = f(u, v)jJ1j. Now discretize SQ into two triangles T0 and T00 , shown in Fig. 1c, then
I¼
ZZ
gðu; v Þ du dv ¼
ZZ T0
SQ
gðu; v Þ du dv þ
ZZ T 00
gðu; v Þ du dv ;
ð5Þ
where T0 = {(u, v):1 6 v 6 1,1 6 u 6 v} and T00 = {(u, v):1 6 u 6 1,u 6 v 6 1}. The second integral of Eq. (5) is then transformed into an integral over the region of the same standard triangle, T0 = {(u, v):1 6 v 6 1, 1 6 u 6 v}, using simple coordinate transformation as shown in Fig. 2, such that:
ZZ
T 00
gðu; v Þ du dv ¼
ZZ
T0
gðv ; uÞ du dv :
ð6Þ
Then Eq. (5) leads us
I¼
ZZ
gðu; v Þ du dv ¼ SQ
ZZ T0
½gðu; v Þ þ gðv ; uÞ du dv ¼
Z
1
1
Z v
Fðu; v Þ du dv ;
ð7Þ
1
where,F(u, v) = g(u, v) + g(v, u). The integral I of Eq. (7) can be further transformed into an integral over the standard 2-square, {(n, g):1 6 n, g 6 1} using standard quadrilateral basis functions, Qi(n, g), as shown in Fig. 3. Assume that
u¼
4 X
ui Q i ðn; gÞ ¼
1 ð1 þ 3n g ngÞ ¼ uðn; gÞ; 4
ð8aÞ
v i Q i ðn; gÞ ¼
1 ð1 n þ 3g ngÞ ¼ v ðn; gÞ; 4
ð8bÞ
i¼1
v¼
4 X i¼1
and
J2 ¼
ou ov ou ov 1 ¼ ð2 n gÞ: on og og on 4
ð8cÞ
Notice that J1 depends on the vertices of the given arbitrary quadrilateral region, but J2 is fixed.
( − 1 ,1 ) 4
( − 1 ,1 ) 4
v
η
3 ( 1, 1 )
3 ( 0,0 )
u
(−1,−1) 1
2 (1,−1)
(a) Standard triangle, ST
ξ
( − 1, − 1) 1
2 (1, − 1)
(b) Standard square, SQ
Fig. 3. Transformation of standard triangle ST into 2-square SQ. (a) Standard triangle, ST. (b) Standard square, SQ.
518
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
y 4 (1 , 4 )
( − 1, 2 ) 1
3 (3,3)
R1
2 (2,1)
x Fig. 4. Quadrilateral region R1.
Then the Eq. (7) reduces to
I¼
Z
1
1
¼
Z
1
Z Z
1
1
Fðuðn; gÞ; v ðn; gÞÞjJ 2 j dn dg 1 1
F 1
1 1 1 ð1 þ 3n g ngÞ; ð1 n þ 3g ngÞ ð2 n gÞ dn dg: 4 4 4
ð9Þ
Now Eq. (9) represents an integral over the standard 2-square region: {(n, g):1 6 n, g 6 1}. Hence using conventional Gauss Legendre quadrature rule for the integral I of Eq. (9), we have
I¼
s X s X 1 1 1 ð2 gj ni Þwi wj F ð1 þ 3ni gj ni gj Þ; ð1 ni þ 3gj ni gj Þ ; 4 4 4 i¼1 j¼1
ð10Þ
where (ni, gj) are Gaussian points in the (n, g) directions of order s, and wi, wj are the corresponding weight coefficients [4]. We can write Eq. (10) as:
I¼
N¼ss X
c0k Fðx0k ; y0k Þ;
ð11aÞ
k¼1
where c0k ; x0k and y0k can be written in the form:
c0k ¼
1 ð2 ni gj Þwi wj ; 4
x0k ¼
1 ð1 þ 3ni gj ni gj Þ; 4
y0k ¼
1 ð1 ni þ 3gj ni gj Þ; 4
ð11bÞ
where k = 1, 2, . . . , N, and i, j = 1, 2, . . . , s The weighting coefficients c0k and sampling points ðx0k ; y0k Þ of various order can be now easily computed from Eq. (11). Using the given program in Mathematica, the outputs of c0k ; x0k and y0k for s = 2, 3, 4, 5, 6,7 are given in Table 1. Now we discretize ST = {(u, v):1 6 v 6 1,1 6 u 6 v} in (u, v)-space of Eq. (7) into 4(n n) = 4n2 right isosceles triangle, each Ti of area 1/(2n2) [13]. Then Eq. (7) reduces to
I¼
4ðnnÞ X Z i¼1
Z
Fðu; v Þ du dv :
ð12Þ
Ti
Since each Ti is to be transformed again into a standard triangle, and using composite integration rule [13] we can obtain the following:
I¼
X 1 N¼ss c0 Hðx0k ; y0k Þ; 2 4n k¼1 k
ð13Þ
where
Hðx0k ; y0k Þ ¼
0 xk þ 2ði nÞ þ 1 y0k þ 2ðj nÞ þ 1 ; 2n 2n i¼0 j¼0 2n2 2n2i 0 0 X X xk þ 2ði nÞ þ 1 yk þ 2ðj nÞ þ 1 F ; ; þ 2n 2n i¼0 j¼0 2n1 X 2n1i X
F
ð14aÞ
and
1 1 ð2 np gq Þwp wq ; x0k ¼ ð1 þ 3np gq np gq Þ; 4 4 1 0 yk ¼ ð1 np þ 3gq np gq Þ; ðk ¼ 1; 2; . . . ; N; p; q ¼ 1; 2; . . . ; sÞ: 4 c0k ¼
ð14bÞ
519
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524 Table 1 Outputs of c0k ; x0k and y0k of Eqs. (11) using the given program in the next section. k
c0k
x0k
y0k
1 2 3 4
Order of Gauss Legendre quadrature rule, s = 2 0.500000000000000 0.211324865405187 0.788675134594813 0.500000000000000
0.744016935856293 0.044658198738520 0.622008467928146 0.410683602522959
0.410683602522959 0.044658198738520 0.622008467928146 0.744016935856293
1 2 3 4 5 6 7 8 9
Order of Gauss Legendre quadrature rule, s = 3 0.154320987654321 0.151284361822039 0.034784464623228 0.342542798671788 0.395061728395062 0.151284361822039 0.273857510685414 0.342542798671788 0.154320987654321
0.874596669241483 0.443649167310371 0.012701665379258 0.830947501931112 0.250000000000000 0.330947501931112 0.787298334620741 0.056350832689629 0.674596669241483
0.674596669241483 0.330947501931112 0.012701665379258 0.056350832689629 0.250000000000000 0.443649167310371 0.787298334620741 0.830947501931112 0.874596669241483
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Order of Gauss Legendre quadrature rule, s = 4 0.060501496642801 0.083869666513367 0.045307001847492 0.008401460977899 0.142982185338485 0.212646651505347 0.140350821011734 0.045307001847492 0.181544850004359 0.284942481998960 0.212646651505347 0.083869666513367 0.112601532307703 0.181544850004359 0.142982185338485 0.060501496642801
0.925747374807601 0.647077355116015 0.283490800681011 0.004820780989426 0.907654989120614 0.561084266085594 0.108906255706834 0.237664467328186 0.884049478270466 0.448887299291690 0.118877821084118 0.554040000062894 0.865957092583479 0.362894210261269 0.293462366058295 0.796525248380506
0.796525248380506 0.554040000062894 0.237664467328186 0.004820780989426 0.293462366058295 0.118877821084118 0.108906255706834 0.283490800681011 0.362894210261269 0.448887299291690 0.561084266085594 0.647077355116015 0.865957092583479 0.884049478270466 0.907654989120614 0.925747374807600
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Order of Gauss Legendre quadrature rule, s = 5 0.028067174431214 0.046275406309135 0.036857657161022 0.015744196426143 0.002633266629203 0.067124593690865 0.114542702111995 0.099488780941957 0.052864972328108 0.015744196426143 0.097927415226499 0.172797751608793 0.161817283950617 0.099488780941957 0.036857657161022 0.097655803573857 0.176220431895882 0.172797751608793 0.114542702111995 0.046275406309135 0.053501082233226 0.097655803573857 0.097927415226499 0.067124593690865 0.028067174431214
0.950889367642309 0.758409434945362 0.476544961484666 0.194680488023970 0.002200555327023 0.942264702861852 0.715982010624261 0.384617327526421 0.053252644428581 0.173030047809011 0.929634884453998 0.653851982579262 0.250000000000000 0.153851982579262 0.429634884453998 0.917005066046144 0.591721954534264 0.115382672473579 0.360956609587106 0.686239721098985 0.908380401265687 0.549294530213163 0.023455038515334 0.502384453182495 0.861470324235019
0.861470324235019 0.686239721098985 0.429634884453998 0.173030047809011 0.002200555327023 0.502384453182495 0.360956609587106 0.153851982579262 0.053252644428581 0.194680488023970 0.023455038515334 0.115382672473579 0.250000000000000 0.384617327526421 0.476544961484666 0.549294530213163 0.591721954534264 0.653851982579262 0.715982010624261 0.758409434945362 0.908380401265687 0.917005066046144 0.929634884453998 0.942264702861852 0.950889367642309
1 2 3 4 5 6 7 8
Order of Gauss Legendre quadrature rule, s = 6 0.014676040844490 0.026712183112376 0.026176910420220 0.016612442896900 0.006278401847429 0.000991080167803 0.035095110260008 0.065074456294084
0.965094665473587 0.824885019554296 0.606455488981548 0.359779268120028 0.141349737547280 0.001140091627990 0.960515083422740 0.801909923278491
0.899844362932718 0.768793881115121 0.564633211304801 0.334071059999927 0.129910390189607 0.001140091627990 0.633163817246677 0.520508849654039 (continued on next page)
520
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
Table 1 (continued) k
c0k
x0k
y0k
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
0.066568810067974 0.046428710417516 0.022046614973247 0.006278401847429 0.053988206897587 0.102236557019614 0.109471725083648 0.083349671145065 0.046428710417516 0.016612442896900 0.063552674420906 0.122376656670072 0.135593779022232 0.109471725083648 0.066568810067974 0.026176910420220 0.055528891524955 0.108102297614921 0.122376656670072 0.102236557019614 0.065074456294084 0.026712183112376 0.028361001521177 0.055528891524955 0.063552674420906 0.053988206897587 0.035095110260008 0.014676040844490
0.554822424771676 0.275782268461456 0.028694769954641 0.129910390189607 0.953380653041526 0.766117524963210 0.474384407091445 0.144925185950153 0.146807931921612 0.334071059999927 0.945323618263202 0.725696554736187 0.383544372033350 0.002853965074949 0.345006147777786 0.564633211304801 0.938189187881988 0.689904156420906 0.303106354353119 0.133711047586252 0.520508849654039 0.768793881115121 0.933609605831142 0.666929060145101 0.251473290143247 0.217708047244823 0.633163817246677 0.899844362932718
0.345006147777786 0.146807931921612 0.028694769954641 0.141349737547280 0.217708047244823 0.133711047586252 0.002853965074949 0.144925185950153 0.275782268461456 0.359779268120028 0.251473290143247 0.303106354353119 0.383544372033350 0.474384407091445 0.554822424771676 0.606455488981548 0.666929060145101 0.689904156420906 0.725696554736187 0.766117524963210 0.801909923278491 0.824885019554296 0.933609605831142 0.938189187881988 0.945323618263202 0.953380653041526 0.960515083422740 0.965094665473587
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
Order of Gauss Legendre quadrature rule, s = 7 0.008383178231877 0.016229336623041 0.018005727931648 0.014218420393005 0.007972981898571 0.002801080689151 0.000426637441423 0.019988306531199 0.039117553014085 0.044437151163993 0.036780461704774 0.022765035551479 0.010110667549803 0.002801080689151 0.031435523240266 0.062362772594912 0.072897093734371 0.063602544468903 0.043312161692773 0.022765035551479 0.007972981898571 0.039901010364923 0.080124975391153 0.095986831742216 0.087344939608496 0.063602544468903 0.036780461704774 0.014218420393005 0.041468269273342 0.084034888207426 0.102482025775969 0.095986831742216 0.072897093734371 0.044437151163993 0.018005727931648 0.033416562465088 0.068124438478366 0.084034888207426 0.080124975391153 0.062362772594912 0.039117553014085 0.016229336623041
0.973906455024851 0.867477088409913 0.695363130529180 0.487276978085690 0.279190825642200 0.107076867761467 0.000647501146528 0.971265451781595 0.854064060795283 0.664529950865235 0.435382796399849 0.206235641934462 0.016701532004414 0.100499858981898 0.966994511011860 0.832372967976233 0.614667579653261 0.351461287844349 0.088254996035437 0.129450392287535 0.264071935323162 0.961830934257069 0.806148389199546 0.554383863533048 0.250000000000000 0.054383863533048 0.306148389199546 0.461830934257069 0.956667357502278 0.779923810422858 0.494100147412834 0.148538712155651 0.197022723101533 0.482846386111557 0.659589933190977 0.952396416732544 0.758232717603808 0.444237776200861 0.064617203600152 0.315003369000558 0.628998310403505 0.823162009532241
0.924309369660667 0.823162009532241 0.659589933190977 0.461830934257069 0.264071935323162 0.100499858981898 0.000647501146528 0.719373646160558 0.628998310403505 0.482846386111557 0.306148389199546 0.129450392287535 0.016701532004414 0.107076867761467 0.387958552708296 0.315003369000558 0.197022723101533 0.054383863533048 0.088254996035437 0.206235641934462 0.279190825642200 0.012723021914310 0.064617203600152 0.148538712155651 0.250000000000000 0.351461287844349 0.435382796399849 0.487276978085690 0.413404596536916 0.444237776200861 0.494100147412834 0.554383863533048 0.614667579653261 0.664529950865235 0.695363130529179 0.744819689989179 0.758232717603808 0.779923810422858 0.806148389199546 0.832372967976233 0.854064060795283 0.867477088409913
521
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524 Table 1 (continued) k
c0k
x0k
y0k
43 44 45 46 47 48 49
0.016339719022330 0.033416562465088 0.041468269273342 0.039901010364923 0.031435523240266 0.019988306531199 0.008383178231877
0.949755413489287 0.744819689989179 0.413404596536916 0.012723021914310 0.387958552708296 0.719373646160558 0.924309369660667
0.949755413489287 0.952396416732544 0.956667357502278 0.961830934257069 0.966994511011860 0.971265451781595 0.973906455024851
3. Algorithm and program Now we give an algorithm and the corresponding Mathematica program of the formulation of the above results in the following: Algorithm. Step 1: Input (conventional) Gauss Legendre quadrature rule of order s. Step 2: Set up a subalgorithm for generating new weighting coefficients c0k and sampling points ðx0k ; y0k Þ to be called STriangle(s) that accepts as input the order of Gauss Legendre quadrature rule s s, and defined by the following: For i = 1, 2, 3, . . . , s For j = 1, 2, 3, . . . , s Set k = i + (s j)s Set the Eq. (11b) Step 3: Call subprogram STriangle(s). Step 4: Input the value of n then T will be discretized into 4 n2 subtriangles Ti. Step 5: Set F(u, v) = f(u, v) + f(v, u). Step 6: For m = 1, 2, 3, . . . , n do steps 7 and 8. Step 7: For k = 1, 2, 3, . . . , s2 Set the Eqs. (14) Step 8: Set the Eq. (13) Step 9: Output for m = 1, 2, 3, . . . , n, the results of the required double integral using 4 m2 subtriangles. Mathematica program < NumericalMath ‘GaussianQuadrature’ STriangle[s0_]:¼Module[{s=s0}, C1 = Table[0, {i, 1, s2}]; X1 = Table[0, {i, 1, s2}]; Y1 = Table[{0, i, 1, s2}]; XI = Table[0, {i, 1, s}]; WI = Table[0, {i, 1, s}]; GQ = GaussianQuadratureWeights [s, 1, 1]; Do [{XIsit = GQsi, 1t; WIsit = GQsi, 2t;}, {i, 1, s}]; For [i = 1, i 6 s, i++, For[j = 1, j 6 s, j++, k = i + (s j) s; 2XI sit XI sjt WI sit WI sjt ; C1 skt ¼ 4 1þ3XI sit XI sjt XI sit XI sjt ; X1 skt ¼ 4 1XI sit þ3XI sjt XI sit XI sjt ; ; ; ; Y1 skt ¼ 4 s = Input[‘‘Enter the Order of Gauss Legendre Quadrature rule:”]; STriangle[s]; TableForm[Table[{PaddedForm[C1sit,{20, 15}], PaddedForm[X1sit//N, {20, 15}], PaddedForm[Y1sit//N, {20, 15}]}, {i, 1, s2}], TableHeadings -> {None,{”ntntnt” ‘‘Ck”,”ntntnt” ‘‘Xk”, ”ntntnt” ‘‘Yk”}}] (*This generates new weighting coefficients and sampling points*) s = Input[”Enter the Order of Gauss Legendre Quadrature rule:”]; STriangle[s] ; c2 ¼ 328125 ; c3 ¼ 239062 ; d1 ¼ 1; d2 ¼ 125 ; d3 ¼ 175 ; c1 ¼ 65625 208 104 208 4 4 8 9 7 6 c1 x þc2 y þc3 x y * * ; ( Example 2 ) f½x ; y :¼ d1 x9 þd2 y7 þd3
F[x_, y_]:¼f[x, y]+f [-y, -x] H1 = Table[0, {i, 1, 4 s2}]; n = Input[‘‘Enter the value of n then discretise T into 4 n2 triangle:”]; App = Table [0, {i, 1, n}];
522
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
Table 2 Evaluation of integrals with present formulation. 4 n2
I1
4 12 4 22 4 32 4 42 4 52 4 62 4 72 4 82 4 92 4 102
Conventional Gauss Legendre quadrature rule, s = 2 3.548632768806294 7.214372270294989 3.549529431843940 14.766944004609170 3.549595019797304 17.750257035191130 3.549607123418959 19.070684778341630 3.549610565443150 19.741189954099270 3.549611827770161 20.116746242388630 3.549612375533267 20.342199117184680 3.549612643456137 20.484773675355270 3.549612786789801 20.578713103462840 3.549612868999844 20.642713389041600
39.039369918293960 42.249315069481930 42.478241857014860 42.523495055881930 42.536817745455030 42.541489797664640 42.543435316001850 42.544391244354330 42.544908789955560 42.545204163254480
4 12 4 22 4 32 4 42 4 52 4 62 4 72 4 82 4 92 4 102
Conventional Gauss Legendre quadrature rule, s = 3 3.549576986770514 14.991195582852550 3.549611923279893 19.267674833364320 3.549612907382378 20.238201270672950 3.549613003510714 20.562029516588160 3.549613020391390 20.693999549518530 3.549613024587089 20.755006700820120 3.549613025900934 20.785863805737210 3.549613026386182 20.802577184929070 3.549613026589042 20.812134862867130 3.549613026682440 20.817848407116390
42.532773862921690 42.542047167997930 42.545962400140290 42.545783237660850 42.545767429494530 42.545763666116570 42.545765986383130 42.545766488461880 42.545766388104710 42.545766345330780
4 12 4 22 4 32 4 42 4 52 4 62 4 72 4 82 4 92 4 102
Conventional Gauss Legendre quadrature rule, s = 4 3.549611375919784 18.644834842838010 3.549613007180398 20.426445648571700 3.549613025681180 20.714072045679120 3.549613026658277 20.787746094199740 3.549613026765570 20.811776994488120 3.549613026783790 20.820951110531610 3.549613026787929 20.824870465865810 3.549613026789089 20.826694184618200 3.549613026789468 20.827602140250910 3.549613026789609 20.828079820693800
42.537308152449840 42.546002670540940 42.545751923101100 42.545764602144810 42.545766495864450 42.545766532601160 42.545766385928470 42.545766357126450 42.545766367740910 42.545766374439800
4 12 4 22 4 32 4 42 4 52 4 62 4 72 4 82 4 92 4 102
Conventional Gauss Legendre quadrature rule, s = 5 3.549612941962180 20.001548496386720 3.549613026386345 20.724196221428660 3.549613026777556 20.806304854241970 3.549613026788828 20.822364616979130 3.549613026789607 20.826571851647100 3.549613026789697 20.827902601602610 3.549613026789713 20.828384024624050 3.549613026789717 20.828576865317250 3.549613026789718 20.828660594976600 3.549613026789717 20.828699421589810
42.543998316928880 42.545752203430210 42.545766387349290 42.545766434803810 42.545766394596410 42.545766356065610 42.545766376187230 42.545766380359400 42.545766378911800 42.545766377750830
4 12 4 22 4 32 4 42 4 52 4 62 4 72 4 82 4 92 4 102
Conventional Gauss Legendre quadrature rule, s = 6 3.549613022130449 20.511614603806180 3.549613026780732 20.801457772983100 3.549613026789571 20.824307659861010 3.549613026789710 20.827736127791060 3.549613026789717 20.828459142086950 3.549613026789717 20.828648914435710 3.549613026789717 20.828707120508940 3.549613026789717 20.828727210885140 3.549613026789717 20.828734824429590 3.549613026789717 20.828737938333210
42.544051506775350 42.545768789991790 42.545766356604670 42.545766313841630 42.545766369364530 42.545766381803220 42.545766378690400 42.545766376692580 42.545766376550800 42.545766376774700
4 12 4 22 4 32 4 42 4 52 4 62 4 72 4 82 4 92 4 102
Conventional Gauss Legendre quadrature rule, s = 7 3.549613026522518 20.707379944169620 3.549613026789507 20.821572098014380 3.549613026789714 20.827856159815450 3.549613026789716 20.828580466438830 3.549613026789717 20.828703557317910 3.549613026789714 20.828730388598700 3.549613026789715 20.828737369893600 3.549613026789716 20.828739447190620 3.549613026789717 20.828740134528660 3.549613026789717 20.828740382547010
42.545475984252790 42.545763813053390 42.545766493516090 42.545766393134410 42.545766374248210 42.545766376430670 42.545766376617350 42.545766376969660 42.545766377156990 42.545766377131950
I2
I3
523
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
For[m = 1, m 6 n, m++, For[k = 1, k 6 s2, k++,
H1
þ 2ðj mÞ þ 1 2m i¼0 j¼0 2m2 2m2i X X X1 skt þ 2ði mÞ þ 1 Y1 skt þ 2ðj mÞ þ 1 þ F ; ; ; 2m 2m i¼0 j¼0 skt
App
¼
smt
2m1 X 2m1i X
¼
X1 F
s2 1 X C1 4m 2 k¼1
skt
þ 2ði mÞ þ 1 Y1 ; 2m
skt
skt H1 skt ; ;
TableForm[Table[{4 i2, PaddedForm[Appsit, {20, 15}]}, {i, 1, n}], TableHeadings -> {None, {‘‘4 n2”, ‘‘ntntnt” ‘‘Result”}}]
4. Numerical examples In this section, we consider three examples to show that the present formulation may be applied to integrate any kind (rational, irrational) of integrals which can not be evaluated even analytically [12]. RR Example 1. Evaluation of I1 ¼ R1 ðx þ yÞ1=2 dy dx, whose exact value is 3.549613026789710, where R1, Shown in Fig. 4, is a quadrilateral region connecting the points (1, 2), (2, 1), (3, 3) and (1, 4). Example 2. Evaluation of I2 ¼
R1 R1 1
1
c1 x8 þc2 y9 þc3 x7 y6 d1 x9 þd2 y7 þd3
dy dx, [12] whose approximate solution is 20.828740382547010 for
; c2 ¼ 328125 ; c3 ¼ 239062 ; d1 ¼ 1; d2 ¼ 125 and d3 ¼ 175 . c1 ¼ 65625 208 104 208 4 4 Example
3. Evaluation
of
42.545766377156990 for c1 ¼
I3 ¼ 65625 ; c2 208
R 1 R 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c1 x16 þ c2 y10 þ c3 x9 y11 dy dx, 1 1 ¼
328125 ; c3 104
¼
[12]
whose
approximate
solution
is
239062 . 208
Now we evaluate the integrals I1, I2 and I3 using the algorithm and the program given in the previous section. The results are shown in Table 2. The same integrals are evaluated also with the conventional Gauss Legendre quadrature rules upto order 30 which are shown in Table 3. Table 3 Evaluation of integrals with conventional Gauss Legendre quadratures. Gauss Legendre quadratures (s)
I1
I2
I3
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
3.538042018537409 3.548681950777742 3.549525572963957 3.549603973886602 3.549612027291594 3.549612911402360 3.549613013031479 3.549613025108831 3.549613026580426 3.549613026763296 3.549613026786368 3.549613026789497 3.549613026789992 3.549613026789129 3.549613026789194 3.549613026786854 3.549613026777713 3.549613026774801 3.549613026767568 3.549613026825007 3.549613026737036 3.549613026569079 3.549613025991078 3.549613031422469 3.549613018504444 3.549613031080794 3.549612941950191 3.549612944461810 3.549612817087430
0.387606580014487 4.041358450537683 10.089910171346000 14.879287405070980 17.674892815377660 19.131995271464990 19.899903477007960 20.322734561053090 20.555472875766140 20.681520680622670 20.749429020941480 20.786044518617900 20.805780742278950 20.816403494092420 20.822114563674550 20.825183168462830 20.826831279051700 20.827716125721150 20.828190998868870 20.828445861309920 20.828582505269080 20.828655837438350 20.828694597062120 20.828718485449080 20.828726175065580 20.828736825240480 20.828658255862420 20.828703859533620 20.828562379650680
14.439843538680000 37.286603723815650 42.573371637903480 42.525315061250620 42.570178907121110 42.550306788523540 42.542625376177950 42.539420680238030 42.542288973666700 42.543761819015370 42.545203032648530 42.545657080679170 42.546038700556750 42.546018388215290 42.545997176492440 42.545891016876810 42.545836220424130 42.545781549438560 42.545761402106100 42.545749236966850 42.545750201699880 42.545753352600670 42.545757858027330 42.545764326372820 42.545761900852400 42.545769038593800 42.545726366048670 42.545734258351530 42.545669325364050
524
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
The results in Table 2 show that the obtained values are stable for large number of triangles (n) with lower order Gauss points (s), while Table 3 shows that if use higher order Gauss quadrature the obtained values oscillate. From Table 3, we may comment that there is no exact order of conventional Gauss quadrature rule those will compute higher order rational and irrational integrals with good accuracy. Both tables also confirm that the present formulation with lower order Gauss points give us better accuracy with high precision compare to the higher order conventional Gauss quadratures. 5. Conclusions In this paper, the formulation of numerical evaluation of double integrals even for rational or irrational expressions over an arbitrary quadrilateral region is discussed in details. At first we have used a coordinate transformation from arbitrary quadrilateral region into a standard 2-square finite element. Then the standard square is discretized into two triangular regions and each of the standard triangle is further discretized into 4 n2 triangles. We then map further each of the standard triangle into the 2-square using standard quadrilateral basis functions. For each triangle we then generate s2 new Gauss points using the lower order conventional Gauss quadrature of order s, and thus the composite numerical integration over the standard triangular finite elements are applied. The program of the present formulations is written by Mathematica, which may be coded easily. We observe that computed results of the given examples converge to the exact solutions correct upto 15 decimal places. We notice that higher order conventional Gauss quadrature rule may produce a larger error than that of a certain lower order, which is shown in Table 3. For this composite numerical integration with lower order Gauss points presented in this paper are more reliable to find the approximate solution that converges to the exact solution. This technique may give fruitful results also for arbitrary triangular regions. Thus the performance of the present formulation is excellent. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]
P.C. Hammer, O.J. Marlowe, A.H. Stroud, Numerical integration over simplexes and cones, Math. Tables Other Aids Comput. 10 (1956) 130–136. P.C. Hammer, A.H. Stroud, Numerical integration over simplexes, Math. Tables Other Aids Comput. 10 (1956) 137–139. P.C. Hammer, A.H. Stroud, Numerical evaluation of multiple integrals, Math. Tables Other Aids Comput. 12 (1958) 272–280. A.H. Stroud, Numerical Quadrature and Solution of Ordinary Differential Equations, Springer-Verlag, New York, Berlin, Heidelberg, 1974. G.R. Cowper, Gaussian quadrature formulas for triangles, Int. J. Numer. Methods Eng. 7 (1973) 405–408. F.G. Lethor, Computation of double integrals over a triangle, J. Comput. Appl. Math. 2 (1976) 219–224. P. Hillion, Numerical integration on a triangle, Int. J. Numer. Methods Eng. 11 (1977) 797–815. M.E. Lauresn, M. Gellert, Some criteria for numerically integrated matrices and quadrature formulas for triangles, Int. J. Numer. Methods Eng. 12 (1978) 67–76. O.C. Zienkiewicz, The Finite Element Method, third ed., McGraw-Hill Inc., 1977. W.B. Bickford, A First Course in the Finite Element Method, Irwin, Illinois, 1990. H.T. Rathod, Some analytical integration formulae for linear isoparametric finite elements, Comput. Struct. 30 (1988) 1101–1109. G. Yagawa, G.W. Ye, S. Yoshimura, A numerical integration scheme for finite element method based on symbolic manipulation, Int. J. Numer. Methods Eng. 29 (1990) 1539–1549. H.T. Rathod, K.V. Nagaraja, B. Venkatesudu, On the application of two symmetric Gauss Legendre quadrature rules for composite numerical integration over a triangular surface, Appl. Math. Comput. 190 (2007) 21–39.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Numerical integrations over an arbitrary quadrilateral region Md. Shafiqul Islam *, M. Alamgir Hossain Department of Mathematics, University of Dhaka, Dhaka 1000, Bangladesh
a r t i c l e
i n f o
Keywords: Double integral Numerical integration Quadrilateral and triangular finite element Gaussian quadrature
a b s t r a c t In this paper, double integrals over an arbitrary quadrilateral are evaluated exploiting finite element method. The physical region is transformed into a standard quadrilateral finite element using the basis functions in local space. Then the standard quadrilateral is subdivided into two triangles, and each triangle is further discretized into 4 n2 right isosceles triangles, with area 2n12 , and thus composite numerical integration is employed. In addition, the affine transformation over each discretized triangle and the use of linearity property of integrals are applied. Finally, each isosceles triangle is transformed into a 2-square finite element to compute new n2 extended symmetric Gauss points and corresponding weight coefficients, where n is the lower order conventional Gauss Legendre quadratures. These new Gauss points and weights are used to compute the double integral. Examples are considered over an arbitrary domain, and rational and irrational integrals which can not be evaluated analytically. Ó 2009 Elsevier Inc. All rights reserved.
1. Introduction Numerical simulation in engineering science and in applied mathematics has become a powerful tool to model the physical phenomena, particularly when analytical solutions are not available and/or are very difficult to obtain. The integrals arising in practical problems are not always simple or polynomial but rational and irrational expressions, in which the quadrature scheme cannot evaluate exactly [12]. Even there is no order of Gauss quadrature that will evaluate these integrals exactly [9,10]. The integration points have to be increased in order to improve the integration accuracy and it is desirable to make these evaluations by using as few Gauss points as possible, from the point of view of the computational efficiency. Among various numerical techniques, the finite element method (FEM) is probably one of the most widely accepted even for the complex geometries. This advantage is supported by the element wise coordinate transformation from one space to the other space. From the literature review we may realize that a lot of works of numerical integration using Gauss quadrature over triangular region has been done [1–8,13], but a limited work is attempted over the quadrilateral region, such as [12]. Rathod [11] presented some analytical formulas for rational integrals over quadrilateral element but it was confined with monomials as numerators. For this, a little work has been done in this study to carry out the development of a good numerical integration technique over the arbitrary convex quadrilateral region with the advent of FEM. Very recently, a rigorous and elaborate survey has been reported in the literature [13] by Rathod et al., and they have derived various orders of extended numerical quadrature rules based on classical Gauss Legendre formula. In their work, a transformation has been used from standard triangular surface to a standard 2-square. All the formulations are derived, and all the examples are tested for certain triangular region. In contrast to this study, we use an arbitrary quadrilateral region with convex and straight sides, which is transformed into a standard square finite element by coordinate transformations. * Corresponding author. E-mail address: mdshafi[email protected] (Md. Shafiqul Islam). 0096-3003/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2009.01.030
516
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
Then the standard square is discretized into two triangular regions and each of these standard triangle is discretized into 4 n2 triangles instead of 3 n2 [13]. We then map further each of the standard triangle into the 2-square using standard quadrilateral basis functions, and taking the sum of the discretized triangles. The subsequent formulations are developed by Mathematica instead of C-Language, which can be coded easily. Examples are considered over an arbitrary domain, and a problem of particular domain which is available in the literature. The results, obtained by the present formulations, converge to the exact solutions correct upto 15 decimal places. 2. Formulation of integrals over an arbitrary quadrilateral region The integral of an arbitrary function, f(x, y) over an arbitrary quadrilateral region AQ is given by
I¼
ZZ
f ðx; yÞ dy dx ¼
ZZ
AQ
f ðx; yÞ dx dy:
ð1Þ
AQ
The integral I of Eq. (1) is then transformed into an integral over the region of the standard quadrilateral, SQ = {(u, v):1 6 u 6 1, 1 6 v 6 1}, shown in Fig. 1, by the changing the coordinates as:
x¼
4 X
xi Q i
and y ¼
i¼1
4 X
yi Q i ;
ð2aÞ
i¼1
where Qi are the bilinear quadrilateral element basis functions in (u, v)-space:
Q 1 ðu; v Þ ¼ ð1 uÞð1 v Þ=4;
Q 2 ðu; v Þ ¼ ð1 þ uÞð1 v Þ=4;
Q 3 ðu; v Þ ¼ ð1 þ uÞð1 þ v Þ=4;
Q 4 ðu; v Þ ¼ ð1 uÞð1 þ v Þ=4:
ð2bÞ
The corresponding Jacobian [11] is then
oðx; yÞ ox oy ox oy ¼ a0 þ a1 u þ a2 v ; ¼ oðu; v Þ ou ov ov ou
ð3aÞ
1 8 1 a1 ¼ ½ðx4 x3 Þðy2 y1 Þ þ ðx1 x2 Þðy4 y3 Þ; 8 1 a2 ¼ ½ðx4 x1 Þðy2 y3 Þ þ ðx3 x2 Þðy4 y1 Þ: 8
ð3bÞ
J1 ¼ where
a0 ¼ ½ðx4 x2 Þðy1 y3 Þ þ ðx3 x1 Þðy4 y2 Þ;
y 3 ( x3 , y3 )
(− 1 , 1 ) 4
v
3 ( 1, 1 )
2 (x2 , y2 )
u
( x4 , y 4 ) 4 1 ( x1 , y1 )
x
( − 1, − 1) 1
(b) Standard square element, SQ
(a) Arbitrary quadrilateral, AQ
v
(− 1 , 1 ) 4
2 (1,−1)
3 ( 1, 1 )
T '' u
T' ( − 1, − 1) 1
2 (1,−1)
(c) Discretize SQ into two triangles T ' and T '' Fig. 1. Transformation of arbitrary quadrilateral (AQ) into standard square (SQ), and discretization SQ into two triangles T0 and T00 .
517
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
v
( − 1 ,1 ) 3
1 ( 1, 1 )
v
( − 1 ,1 ) 3
T '' u
u T
2 (1, − 1)
'
2 (1, − 1)
( − 1, − 1) 1
Fig. 2. Transformation of T00 into T0 .
Now use Eqs. (1)–(3) to obtain,
I¼
ZZ
f ðx; yÞ dy dx ¼
ZZ
AQ
SQ
f ðxðu; v Þ; yðu; v ÞÞjJ 1 j du dv ¼
ZZ
gðu; v Þ du dv ¼ SQ
Z
1 1
Z
1
gðu; v Þ du dv ;
ð4Þ
1
where, g(u, v) = f(u, v)jJ1j. Now discretize SQ into two triangles T0 and T00 , shown in Fig. 1c, then
I¼
ZZ
gðu; v Þ du dv ¼
ZZ T0
SQ
gðu; v Þ du dv þ
ZZ T 00
gðu; v Þ du dv ;
ð5Þ
where T0 = {(u, v):1 6 v 6 1,1 6 u 6 v} and T00 = {(u, v):1 6 u 6 1,u 6 v 6 1}. The second integral of Eq. (5) is then transformed into an integral over the region of the same standard triangle, T0 = {(u, v):1 6 v 6 1, 1 6 u 6 v}, using simple coordinate transformation as shown in Fig. 2, such that:
ZZ
T 00
gðu; v Þ du dv ¼
ZZ
T0
gðv ; uÞ du dv :
ð6Þ
Then Eq. (5) leads us
I¼
ZZ
gðu; v Þ du dv ¼ SQ
ZZ T0
½gðu; v Þ þ gðv ; uÞ du dv ¼
Z
1
1
Z v
Fðu; v Þ du dv ;
ð7Þ
1
where,F(u, v) = g(u, v) + g(v, u). The integral I of Eq. (7) can be further transformed into an integral over the standard 2-square, {(n, g):1 6 n, g 6 1} using standard quadrilateral basis functions, Qi(n, g), as shown in Fig. 3. Assume that
u¼
4 X
ui Q i ðn; gÞ ¼
1 ð1 þ 3n g ngÞ ¼ uðn; gÞ; 4
ð8aÞ
v i Q i ðn; gÞ ¼
1 ð1 n þ 3g ngÞ ¼ v ðn; gÞ; 4
ð8bÞ
i¼1
v¼
4 X i¼1
and
J2 ¼
ou ov ou ov 1 ¼ ð2 n gÞ: on og og on 4
ð8cÞ
Notice that J1 depends on the vertices of the given arbitrary quadrilateral region, but J2 is fixed.
( − 1 ,1 ) 4
( − 1 ,1 ) 4
v
η
3 ( 1, 1 )
3 ( 0,0 )
u
(−1,−1) 1
2 (1,−1)
(a) Standard triangle, ST
ξ
( − 1, − 1) 1
2 (1, − 1)
(b) Standard square, SQ
Fig. 3. Transformation of standard triangle ST into 2-square SQ. (a) Standard triangle, ST. (b) Standard square, SQ.
518
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
y 4 (1 , 4 )
( − 1, 2 ) 1
3 (3,3)
R1
2 (2,1)
x Fig. 4. Quadrilateral region R1.
Then the Eq. (7) reduces to
I¼
Z
1
1
¼
Z
1
Z Z
1
1
Fðuðn; gÞ; v ðn; gÞÞjJ 2 j dn dg 1 1
F 1
1 1 1 ð1 þ 3n g ngÞ; ð1 n þ 3g ngÞ ð2 n gÞ dn dg: 4 4 4
ð9Þ
Now Eq. (9) represents an integral over the standard 2-square region: {(n, g):1 6 n, g 6 1}. Hence using conventional Gauss Legendre quadrature rule for the integral I of Eq. (9), we have
I¼
s X s X 1 1 1 ð2 gj ni Þwi wj F ð1 þ 3ni gj ni gj Þ; ð1 ni þ 3gj ni gj Þ ; 4 4 4 i¼1 j¼1
ð10Þ
where (ni, gj) are Gaussian points in the (n, g) directions of order s, and wi, wj are the corresponding weight coefficients [4]. We can write Eq. (10) as:
I¼
N¼ss X
c0k Fðx0k ; y0k Þ;
ð11aÞ
k¼1
where c0k ; x0k and y0k can be written in the form:
c0k ¼
1 ð2 ni gj Þwi wj ; 4
x0k ¼
1 ð1 þ 3ni gj ni gj Þ; 4
y0k ¼
1 ð1 ni þ 3gj ni gj Þ; 4
ð11bÞ
where k = 1, 2, . . . , N, and i, j = 1, 2, . . . , s The weighting coefficients c0k and sampling points ðx0k ; y0k Þ of various order can be now easily computed from Eq. (11). Using the given program in Mathematica, the outputs of c0k ; x0k and y0k for s = 2, 3, 4, 5, 6,7 are given in Table 1. Now we discretize ST = {(u, v):1 6 v 6 1,1 6 u 6 v} in (u, v)-space of Eq. (7) into 4(n n) = 4n2 right isosceles triangle, each Ti of area 1/(2n2) [13]. Then Eq. (7) reduces to
I¼
4ðnnÞ X Z i¼1
Z
Fðu; v Þ du dv :
ð12Þ
Ti
Since each Ti is to be transformed again into a standard triangle, and using composite integration rule [13] we can obtain the following:
I¼
X 1 N¼ss c0 Hðx0k ; y0k Þ; 2 4n k¼1 k
ð13Þ
where
Hðx0k ; y0k Þ ¼
0 xk þ 2ði nÞ þ 1 y0k þ 2ðj nÞ þ 1 ; 2n 2n i¼0 j¼0 2n2 2n2i 0 0 X X xk þ 2ði nÞ þ 1 yk þ 2ðj nÞ þ 1 F ; ; þ 2n 2n i¼0 j¼0 2n1 X 2n1i X
F
ð14aÞ
and
1 1 ð2 np gq Þwp wq ; x0k ¼ ð1 þ 3np gq np gq Þ; 4 4 1 0 yk ¼ ð1 np þ 3gq np gq Þ; ðk ¼ 1; 2; . . . ; N; p; q ¼ 1; 2; . . . ; sÞ: 4 c0k ¼
ð14bÞ
519
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524 Table 1 Outputs of c0k ; x0k and y0k of Eqs. (11) using the given program in the next section. k
c0k
x0k
y0k
1 2 3 4
Order of Gauss Legendre quadrature rule, s = 2 0.500000000000000 0.211324865405187 0.788675134594813 0.500000000000000
0.744016935856293 0.044658198738520 0.622008467928146 0.410683602522959
0.410683602522959 0.044658198738520 0.622008467928146 0.744016935856293
1 2 3 4 5 6 7 8 9
Order of Gauss Legendre quadrature rule, s = 3 0.154320987654321 0.151284361822039 0.034784464623228 0.342542798671788 0.395061728395062 0.151284361822039 0.273857510685414 0.342542798671788 0.154320987654321
0.874596669241483 0.443649167310371 0.012701665379258 0.830947501931112 0.250000000000000 0.330947501931112 0.787298334620741 0.056350832689629 0.674596669241483
0.674596669241483 0.330947501931112 0.012701665379258 0.056350832689629 0.250000000000000 0.443649167310371 0.787298334620741 0.830947501931112 0.874596669241483
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Order of Gauss Legendre quadrature rule, s = 4 0.060501496642801 0.083869666513367 0.045307001847492 0.008401460977899 0.142982185338485 0.212646651505347 0.140350821011734 0.045307001847492 0.181544850004359 0.284942481998960 0.212646651505347 0.083869666513367 0.112601532307703 0.181544850004359 0.142982185338485 0.060501496642801
0.925747374807601 0.647077355116015 0.283490800681011 0.004820780989426 0.907654989120614 0.561084266085594 0.108906255706834 0.237664467328186 0.884049478270466 0.448887299291690 0.118877821084118 0.554040000062894 0.865957092583479 0.362894210261269 0.293462366058295 0.796525248380506
0.796525248380506 0.554040000062894 0.237664467328186 0.004820780989426 0.293462366058295 0.118877821084118 0.108906255706834 0.283490800681011 0.362894210261269 0.448887299291690 0.561084266085594 0.647077355116015 0.865957092583479 0.884049478270466 0.907654989120614 0.925747374807600
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Order of Gauss Legendre quadrature rule, s = 5 0.028067174431214 0.046275406309135 0.036857657161022 0.015744196426143 0.002633266629203 0.067124593690865 0.114542702111995 0.099488780941957 0.052864972328108 0.015744196426143 0.097927415226499 0.172797751608793 0.161817283950617 0.099488780941957 0.036857657161022 0.097655803573857 0.176220431895882 0.172797751608793 0.114542702111995 0.046275406309135 0.053501082233226 0.097655803573857 0.097927415226499 0.067124593690865 0.028067174431214
0.950889367642309 0.758409434945362 0.476544961484666 0.194680488023970 0.002200555327023 0.942264702861852 0.715982010624261 0.384617327526421 0.053252644428581 0.173030047809011 0.929634884453998 0.653851982579262 0.250000000000000 0.153851982579262 0.429634884453998 0.917005066046144 0.591721954534264 0.115382672473579 0.360956609587106 0.686239721098985 0.908380401265687 0.549294530213163 0.023455038515334 0.502384453182495 0.861470324235019
0.861470324235019 0.686239721098985 0.429634884453998 0.173030047809011 0.002200555327023 0.502384453182495 0.360956609587106 0.153851982579262 0.053252644428581 0.194680488023970 0.023455038515334 0.115382672473579 0.250000000000000 0.384617327526421 0.476544961484666 0.549294530213163 0.591721954534264 0.653851982579262 0.715982010624261 0.758409434945362 0.908380401265687 0.917005066046144 0.929634884453998 0.942264702861852 0.950889367642309
1 2 3 4 5 6 7 8
Order of Gauss Legendre quadrature rule, s = 6 0.014676040844490 0.026712183112376 0.026176910420220 0.016612442896900 0.006278401847429 0.000991080167803 0.035095110260008 0.065074456294084
0.965094665473587 0.824885019554296 0.606455488981548 0.359779268120028 0.141349737547280 0.001140091627990 0.960515083422740 0.801909923278491
0.899844362932718 0.768793881115121 0.564633211304801 0.334071059999927 0.129910390189607 0.001140091627990 0.633163817246677 0.520508849654039 (continued on next page)
520
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
Table 1 (continued) k
c0k
x0k
y0k
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
0.066568810067974 0.046428710417516 0.022046614973247 0.006278401847429 0.053988206897587 0.102236557019614 0.109471725083648 0.083349671145065 0.046428710417516 0.016612442896900 0.063552674420906 0.122376656670072 0.135593779022232 0.109471725083648 0.066568810067974 0.026176910420220 0.055528891524955 0.108102297614921 0.122376656670072 0.102236557019614 0.065074456294084 0.026712183112376 0.028361001521177 0.055528891524955 0.063552674420906 0.053988206897587 0.035095110260008 0.014676040844490
0.554822424771676 0.275782268461456 0.028694769954641 0.129910390189607 0.953380653041526 0.766117524963210 0.474384407091445 0.144925185950153 0.146807931921612 0.334071059999927 0.945323618263202 0.725696554736187 0.383544372033350 0.002853965074949 0.345006147777786 0.564633211304801 0.938189187881988 0.689904156420906 0.303106354353119 0.133711047586252 0.520508849654039 0.768793881115121 0.933609605831142 0.666929060145101 0.251473290143247 0.217708047244823 0.633163817246677 0.899844362932718
0.345006147777786 0.146807931921612 0.028694769954641 0.141349737547280 0.217708047244823 0.133711047586252 0.002853965074949 0.144925185950153 0.275782268461456 0.359779268120028 0.251473290143247 0.303106354353119 0.383544372033350 0.474384407091445 0.554822424771676 0.606455488981548 0.666929060145101 0.689904156420906 0.725696554736187 0.766117524963210 0.801909923278491 0.824885019554296 0.933609605831142 0.938189187881988 0.945323618263202 0.953380653041526 0.960515083422740 0.965094665473587
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
Order of Gauss Legendre quadrature rule, s = 7 0.008383178231877 0.016229336623041 0.018005727931648 0.014218420393005 0.007972981898571 0.002801080689151 0.000426637441423 0.019988306531199 0.039117553014085 0.044437151163993 0.036780461704774 0.022765035551479 0.010110667549803 0.002801080689151 0.031435523240266 0.062362772594912 0.072897093734371 0.063602544468903 0.043312161692773 0.022765035551479 0.007972981898571 0.039901010364923 0.080124975391153 0.095986831742216 0.087344939608496 0.063602544468903 0.036780461704774 0.014218420393005 0.041468269273342 0.084034888207426 0.102482025775969 0.095986831742216 0.072897093734371 0.044437151163993 0.018005727931648 0.033416562465088 0.068124438478366 0.084034888207426 0.080124975391153 0.062362772594912 0.039117553014085 0.016229336623041
0.973906455024851 0.867477088409913 0.695363130529180 0.487276978085690 0.279190825642200 0.107076867761467 0.000647501146528 0.971265451781595 0.854064060795283 0.664529950865235 0.435382796399849 0.206235641934462 0.016701532004414 0.100499858981898 0.966994511011860 0.832372967976233 0.614667579653261 0.351461287844349 0.088254996035437 0.129450392287535 0.264071935323162 0.961830934257069 0.806148389199546 0.554383863533048 0.250000000000000 0.054383863533048 0.306148389199546 0.461830934257069 0.956667357502278 0.779923810422858 0.494100147412834 0.148538712155651 0.197022723101533 0.482846386111557 0.659589933190977 0.952396416732544 0.758232717603808 0.444237776200861 0.064617203600152 0.315003369000558 0.628998310403505 0.823162009532241
0.924309369660667 0.823162009532241 0.659589933190977 0.461830934257069 0.264071935323162 0.100499858981898 0.000647501146528 0.719373646160558 0.628998310403505 0.482846386111557 0.306148389199546 0.129450392287535 0.016701532004414 0.107076867761467 0.387958552708296 0.315003369000558 0.197022723101533 0.054383863533048 0.088254996035437 0.206235641934462 0.279190825642200 0.012723021914310 0.064617203600152 0.148538712155651 0.250000000000000 0.351461287844349 0.435382796399849 0.487276978085690 0.413404596536916 0.444237776200861 0.494100147412834 0.554383863533048 0.614667579653261 0.664529950865235 0.695363130529179 0.744819689989179 0.758232717603808 0.779923810422858 0.806148389199546 0.832372967976233 0.854064060795283 0.867477088409913
521
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524 Table 1 (continued) k
c0k
x0k
y0k
43 44 45 46 47 48 49
0.016339719022330 0.033416562465088 0.041468269273342 0.039901010364923 0.031435523240266 0.019988306531199 0.008383178231877
0.949755413489287 0.744819689989179 0.413404596536916 0.012723021914310 0.387958552708296 0.719373646160558 0.924309369660667
0.949755413489287 0.952396416732544 0.956667357502278 0.961830934257069 0.966994511011860 0.971265451781595 0.973906455024851
3. Algorithm and program Now we give an algorithm and the corresponding Mathematica program of the formulation of the above results in the following: Algorithm. Step 1: Input (conventional) Gauss Legendre quadrature rule of order s. Step 2: Set up a subalgorithm for generating new weighting coefficients c0k and sampling points ðx0k ; y0k Þ to be called STriangle(s) that accepts as input the order of Gauss Legendre quadrature rule s s, and defined by the following: For i = 1, 2, 3, . . . , s For j = 1, 2, 3, . . . , s Set k = i + (s j)s Set the Eq. (11b) Step 3: Call subprogram STriangle(s). Step 4: Input the value of n then T will be discretized into 4 n2 subtriangles Ti. Step 5: Set F(u, v) = f(u, v) + f(v, u). Step 6: For m = 1, 2, 3, . . . , n do steps 7 and 8. Step 7: For k = 1, 2, 3, . . . , s2 Set the Eqs. (14) Step 8: Set the Eq. (13) Step 9: Output for m = 1, 2, 3, . . . , n, the results of the required double integral using 4 m2 subtriangles. Mathematica program < NumericalMath ‘GaussianQuadrature’ STriangle[s0_]:¼Module[{s=s0}, C1 = Table[0, {i, 1, s2}]; X1 = Table[0, {i, 1, s2}]; Y1 = Table[{0, i, 1, s2}]; XI = Table[0, {i, 1, s}]; WI = Table[0, {i, 1, s}]; GQ = GaussianQuadratureWeights [s, 1, 1]; Do [{XIsit = GQsi, 1t; WIsit = GQsi, 2t;}, {i, 1, s}]; For [i = 1, i 6 s, i++, For[j = 1, j 6 s, j++, k = i + (s j) s; 2XI sit XI sjt WI sit WI sjt ; C1 skt ¼ 4 1þ3XI sit XI sjt XI sit XI sjt ; X1 skt ¼ 4 1XI sit þ3XI sjt XI sit XI sjt ; ; ; ; Y1 skt ¼ 4 s = Input[‘‘Enter the Order of Gauss Legendre Quadrature rule:”]; STriangle[s]; TableForm[Table[{PaddedForm[C1sit,{20, 15}], PaddedForm[X1sit//N, {20, 15}], PaddedForm[Y1sit//N, {20, 15}]}, {i, 1, s2}], TableHeadings -> {None,{”ntntnt” ‘‘Ck”,”ntntnt” ‘‘Xk”, ”ntntnt” ‘‘Yk”}}] (*This generates new weighting coefficients and sampling points*) s = Input[”Enter the Order of Gauss Legendre Quadrature rule:”]; STriangle[s] ; c2 ¼ 328125 ; c3 ¼ 239062 ; d1 ¼ 1; d2 ¼ 125 ; d3 ¼ 175 ; c1 ¼ 65625 208 104 208 4 4 8 9 7 6 c1 x þc2 y þc3 x y * * ; ( Example 2 ) f½x ; y :¼ d1 x9 þd2 y7 þd3
F[x_, y_]:¼f[x, y]+f [-y, -x] H1 = Table[0, {i, 1, 4 s2}]; n = Input[‘‘Enter the value of n then discretise T into 4 n2 triangle:”]; App = Table [0, {i, 1, n}];
522
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
Table 2 Evaluation of integrals with present formulation. 4 n2
I1
4 12 4 22 4 32 4 42 4 52 4 62 4 72 4 82 4 92 4 102
Conventional Gauss Legendre quadrature rule, s = 2 3.548632768806294 7.214372270294989 3.549529431843940 14.766944004609170 3.549595019797304 17.750257035191130 3.549607123418959 19.070684778341630 3.549610565443150 19.741189954099270 3.549611827770161 20.116746242388630 3.549612375533267 20.342199117184680 3.549612643456137 20.484773675355270 3.549612786789801 20.578713103462840 3.549612868999844 20.642713389041600
39.039369918293960 42.249315069481930 42.478241857014860 42.523495055881930 42.536817745455030 42.541489797664640 42.543435316001850 42.544391244354330 42.544908789955560 42.545204163254480
4 12 4 22 4 32 4 42 4 52 4 62 4 72 4 82 4 92 4 102
Conventional Gauss Legendre quadrature rule, s = 3 3.549576986770514 14.991195582852550 3.549611923279893 19.267674833364320 3.549612907382378 20.238201270672950 3.549613003510714 20.562029516588160 3.549613020391390 20.693999549518530 3.549613024587089 20.755006700820120 3.549613025900934 20.785863805737210 3.549613026386182 20.802577184929070 3.549613026589042 20.812134862867130 3.549613026682440 20.817848407116390
42.532773862921690 42.542047167997930 42.545962400140290 42.545783237660850 42.545767429494530 42.545763666116570 42.545765986383130 42.545766488461880 42.545766388104710 42.545766345330780
4 12 4 22 4 32 4 42 4 52 4 62 4 72 4 82 4 92 4 102
Conventional Gauss Legendre quadrature rule, s = 4 3.549611375919784 18.644834842838010 3.549613007180398 20.426445648571700 3.549613025681180 20.714072045679120 3.549613026658277 20.787746094199740 3.549613026765570 20.811776994488120 3.549613026783790 20.820951110531610 3.549613026787929 20.824870465865810 3.549613026789089 20.826694184618200 3.549613026789468 20.827602140250910 3.549613026789609 20.828079820693800
42.537308152449840 42.546002670540940 42.545751923101100 42.545764602144810 42.545766495864450 42.545766532601160 42.545766385928470 42.545766357126450 42.545766367740910 42.545766374439800
4 12 4 22 4 32 4 42 4 52 4 62 4 72 4 82 4 92 4 102
Conventional Gauss Legendre quadrature rule, s = 5 3.549612941962180 20.001548496386720 3.549613026386345 20.724196221428660 3.549613026777556 20.806304854241970 3.549613026788828 20.822364616979130 3.549613026789607 20.826571851647100 3.549613026789697 20.827902601602610 3.549613026789713 20.828384024624050 3.549613026789717 20.828576865317250 3.549613026789718 20.828660594976600 3.549613026789717 20.828699421589810
42.543998316928880 42.545752203430210 42.545766387349290 42.545766434803810 42.545766394596410 42.545766356065610 42.545766376187230 42.545766380359400 42.545766378911800 42.545766377750830
4 12 4 22 4 32 4 42 4 52 4 62 4 72 4 82 4 92 4 102
Conventional Gauss Legendre quadrature rule, s = 6 3.549613022130449 20.511614603806180 3.549613026780732 20.801457772983100 3.549613026789571 20.824307659861010 3.549613026789710 20.827736127791060 3.549613026789717 20.828459142086950 3.549613026789717 20.828648914435710 3.549613026789717 20.828707120508940 3.549613026789717 20.828727210885140 3.549613026789717 20.828734824429590 3.549613026789717 20.828737938333210
42.544051506775350 42.545768789991790 42.545766356604670 42.545766313841630 42.545766369364530 42.545766381803220 42.545766378690400 42.545766376692580 42.545766376550800 42.545766376774700
4 12 4 22 4 32 4 42 4 52 4 62 4 72 4 82 4 92 4 102
Conventional Gauss Legendre quadrature rule, s = 7 3.549613026522518 20.707379944169620 3.549613026789507 20.821572098014380 3.549613026789714 20.827856159815450 3.549613026789716 20.828580466438830 3.549613026789717 20.828703557317910 3.549613026789714 20.828730388598700 3.549613026789715 20.828737369893600 3.549613026789716 20.828739447190620 3.549613026789717 20.828740134528660 3.549613026789717 20.828740382547010
42.545475984252790 42.545763813053390 42.545766493516090 42.545766393134410 42.545766374248210 42.545766376430670 42.545766376617350 42.545766376969660 42.545766377156990 42.545766377131950
I2
I3
523
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
For[m = 1, m 6 n, m++, For[k = 1, k 6 s2, k++,
H1
þ 2ðj mÞ þ 1 2m i¼0 j¼0 2m2 2m2i X X X1 skt þ 2ði mÞ þ 1 Y1 skt þ 2ðj mÞ þ 1 þ F ; ; ; 2m 2m i¼0 j¼0 skt
App
¼
smt
2m1 X 2m1i X
¼
X1 F
s2 1 X C1 4m 2 k¼1
skt
þ 2ði mÞ þ 1 Y1 ; 2m
skt
skt H1 skt ; ;
TableForm[Table[{4 i2, PaddedForm[Appsit, {20, 15}]}, {i, 1, n}], TableHeadings -> {None, {‘‘4 n2”, ‘‘ntntnt” ‘‘Result”}}]
4. Numerical examples In this section, we consider three examples to show that the present formulation may be applied to integrate any kind (rational, irrational) of integrals which can not be evaluated even analytically [12]. RR Example 1. Evaluation of I1 ¼ R1 ðx þ yÞ1=2 dy dx, whose exact value is 3.549613026789710, where R1, Shown in Fig. 4, is a quadrilateral region connecting the points (1, 2), (2, 1), (3, 3) and (1, 4). Example 2. Evaluation of I2 ¼
R1 R1 1
1
c1 x8 þc2 y9 þc3 x7 y6 d1 x9 þd2 y7 þd3
dy dx, [12] whose approximate solution is 20.828740382547010 for
; c2 ¼ 328125 ; c3 ¼ 239062 ; d1 ¼ 1; d2 ¼ 125 and d3 ¼ 175 . c1 ¼ 65625 208 104 208 4 4 Example
3. Evaluation
of
42.545766377156990 for c1 ¼
I3 ¼ 65625 ; c2 208
R 1 R 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c1 x16 þ c2 y10 þ c3 x9 y11 dy dx, 1 1 ¼
328125 ; c3 104
¼
[12]
whose
approximate
solution
is
239062 . 208
Now we evaluate the integrals I1, I2 and I3 using the algorithm and the program given in the previous section. The results are shown in Table 2. The same integrals are evaluated also with the conventional Gauss Legendre quadrature rules upto order 30 which are shown in Table 3. Table 3 Evaluation of integrals with conventional Gauss Legendre quadratures. Gauss Legendre quadratures (s)
I1
I2
I3
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
3.538042018537409 3.548681950777742 3.549525572963957 3.549603973886602 3.549612027291594 3.549612911402360 3.549613013031479 3.549613025108831 3.549613026580426 3.549613026763296 3.549613026786368 3.549613026789497 3.549613026789992 3.549613026789129 3.549613026789194 3.549613026786854 3.549613026777713 3.549613026774801 3.549613026767568 3.549613026825007 3.549613026737036 3.549613026569079 3.549613025991078 3.549613031422469 3.549613018504444 3.549613031080794 3.549612941950191 3.549612944461810 3.549612817087430
0.387606580014487 4.041358450537683 10.089910171346000 14.879287405070980 17.674892815377660 19.131995271464990 19.899903477007960 20.322734561053090 20.555472875766140 20.681520680622670 20.749429020941480 20.786044518617900 20.805780742278950 20.816403494092420 20.822114563674550 20.825183168462830 20.826831279051700 20.827716125721150 20.828190998868870 20.828445861309920 20.828582505269080 20.828655837438350 20.828694597062120 20.828718485449080 20.828726175065580 20.828736825240480 20.828658255862420 20.828703859533620 20.828562379650680
14.439843538680000 37.286603723815650 42.573371637903480 42.525315061250620 42.570178907121110 42.550306788523540 42.542625376177950 42.539420680238030 42.542288973666700 42.543761819015370 42.545203032648530 42.545657080679170 42.546038700556750 42.546018388215290 42.545997176492440 42.545891016876810 42.545836220424130 42.545781549438560 42.545761402106100 42.545749236966850 42.545750201699880 42.545753352600670 42.545757858027330 42.545764326372820 42.545761900852400 42.545769038593800 42.545726366048670 42.545734258351530 42.545669325364050
524
Md. Shafiqul Islam, M. Alamgir Hossain / Applied Mathematics and Computation 210 (2009) 515–524
The results in Table 2 show that the obtained values are stable for large number of triangles (n) with lower order Gauss points (s), while Table 3 shows that if use higher order Gauss quadrature the obtained values oscillate. From Table 3, we may comment that there is no exact order of conventional Gauss quadrature rule those will compute higher order rational and irrational integrals with good accuracy. Both tables also confirm that the present formulation with lower order Gauss points give us better accuracy with high precision compare to the higher order conventional Gauss quadratures. 5. Conclusions In this paper, the formulation of numerical evaluation of double integrals even for rational or irrational expressions over an arbitrary quadrilateral region is discussed in details. At first we have used a coordinate transformation from arbitrary quadrilateral region into a standard 2-square finite element. Then the standard square is discretized into two triangular regions and each of the standard triangle is further discretized into 4 n2 triangles. We then map further each of the standard triangle into the 2-square using standard quadrilateral basis functions. For each triangle we then generate s2 new Gauss points using the lower order conventional Gauss quadrature of order s, and thus the composite numerical integration over the standard triangular finite elements are applied. The program of the present formulations is written by Mathematica, which may be coded easily. We observe that computed results of the given examples converge to the exact solutions correct upto 15 decimal places. We notice that higher order conventional Gauss quadrature rule may produce a larger error than that of a certain lower order, which is shown in Table 3. For this composite numerical integration with lower order Gauss points presented in this paper are more reliable to find the approximate solution that converges to the exact solution. This technique may give fruitful results also for arbitrary triangular regions. Thus the performance of the present formulation is excellent. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]
P.C. Hammer, O.J. Marlowe, A.H. Stroud, Numerical integration over simplexes and cones, Math. Tables Other Aids Comput. 10 (1956) 130–136. P.C. Hammer, A.H. Stroud, Numerical integration over simplexes, Math. Tables Other Aids Comput. 10 (1956) 137–139. P.C. Hammer, A.H. Stroud, Numerical evaluation of multiple integrals, Math. Tables Other Aids Comput. 12 (1958) 272–280. A.H. Stroud, Numerical Quadrature and Solution of Ordinary Differential Equations, Springer-Verlag, New York, Berlin, Heidelberg, 1974. G.R. Cowper, Gaussian quadrature formulas for triangles, Int. J. Numer. Methods Eng. 7 (1973) 405–408. F.G. Lethor, Computation of double integrals over a triangle, J. Comput. Appl. Math. 2 (1976) 219–224. P. Hillion, Numerical integration on a triangle, Int. J. Numer. Methods Eng. 11 (1977) 797–815. M.E. Lauresn, M. Gellert, Some criteria for numerically integrated matrices and quadrature formulas for triangles, Int. J. Numer. Methods Eng. 12 (1978) 67–76. O.C. Zienkiewicz, The Finite Element Method, third ed., McGraw-Hill Inc., 1977. W.B. Bickford, A First Course in the Finite Element Method, Irwin, Illinois, 1990. H.T. Rathod, Some analytical integration formulae for linear isoparametric finite elements, Comput. Struct. 30 (1988) 1101–1109. G. Yagawa, G.W. Ye, S. Yoshimura, A numerical integration scheme for finite element method based on symbolic manipulation, Int. J. Numer. Methods Eng. 29 (1990) 1539–1549. H.T. Rathod, K.V. Nagaraja, B. Venkatesudu, On the application of two symmetric Gauss Legendre quadrature rules for composite numerical integration over a triangular surface, Appl. Math. Comput. 190 (2007) 21–39.
