LIBRARY OF THE 
 
 UNIVERSITY OF ILLINOIS 
 
 AT URBANA-CHAMPAIGN 
 
 510.84 
 I46r 
 
 no. 93-98 
 cop 3 
 
Digitized by the Internet Archive 
 
 in 2013 
 
 http://archive.org/details/numericalquadrat95mill 
 
5JO-84- 
 
 rvo.95" 
 
 5 
 
 cop. 
 
 UNIVERSITY OF ILLINOIS 
 GRADUATE COLLEGE 
 DIGITAL COMPUTER LABORATORY 
 
 REPORT NO. 95 
 
 NUMERICAL QUADRATURE OVER A RECTANGULAR DOMAIN 
 IN TWO OR MORE DIMENSIONS 
 
 by 
 
 J. Co P. Miller 
 December 11, 1959 
 
 (This is the third of a series of papers, all of which 
 will be published in Mathematics of Computation) 
 
 This work was supported in part by the Office of 
 Naval Research under Contract Nonr- 183^(27) « 
 
Numerical Qua.drature over a, Rectangular Domain 
 in Two or More Dimensions 
 
 By J.C.P. Miller 
 III Quadrature of a. Harmonic Integrand 
 
 Introductory 
 
 1. In Note 1(1), §5, formula (B),<^7, formula (B'), and in § 9, also in Note Il(l) 
 in several places, we have seen how the error term is very much reduced if the 
 integrand f(x,y) is a harmonic function, that is, if w f - 0. In this note we 
 pursue further this special case, in which especially high accuracy is attainable 
 with few points. 
 
 It may not be often that the integrand will have this special form, but it 
 seems worthwhile to develop a few of the interesting formulae. We start by 
 obtaining expansions for n variables, a.nd more extensive ones for two variables, 
 and then obta.in and consider special quadrature formulae. 
 
 2. Expansions . As in 11(1), §2, ve develop f(^, V ...,xJ as a. Taylor series 
 in even powers of ea.ch of the variables Xj,. Then, using v f n = whenever it 
 
 is applicable, we obtain 
 
 r j = l/(2h) n = (2h)~ n /j^ J f( Xl ,x 2 ,...,X n ) dx 1 , 
 
 dx_ . . . dx 
 2 n 
 
 . f + ^ itj3 f + *L i6^6 + h! ( l6 8 192^ 
 o + 5.* 3^0?! 3 J 9.' v 5 ** 5 <*- ' o 
 
 h 10 ,128 _> 6 1280 10. 
 
 (2.1^ + ttt(— rf a +± *t<P )f o 
 
 2 r 64 12 U096 12 61181^ 4 ^8 707384 /L2. 
 
 k 12 
 + 13. 
 
 V 
 
 where, as before, extended, 
 
 6 8 
 
 (2:2) nr\ - Z jV,, Cr\ - s / f ° a Jl\ - z g g /° 2 g 
 
 » X r 3x s ^ r ^ x t 8x r 8x s 3 \ * x u 
 
 etc., the summations extending over all possible combinations of r, s, t, ...with 
 no two equal. 
 (l) J.C.P.Miller,M.T.A.C, in the press. 
 
- 2 - 
 
 We have likewise, labelling the symmetrical sets of points as in Note II, 
 the expansions for sums of values of f over the sets 
 
 (2.31) f Q 
 
 (2.32) a(a) 2nf Q - ***- tft Q + *£- rft Q + ^(^^ 
 
 -10h 10 a 10 '■ - 1212 
 
 10 
 
 B»V/V - ^W^c^-e^/V--- 
 
 (2.33) P (b) 2n(n-l)f -^Mn-10K?f - i 2 ^ U-i6) J 6 f 
 
 IP IP 
 
 > ^, f 2(n-34)^ 12 -3(n + 362)C7 L2 -6(n-728)^ < ^-f 6.(n-102^ 12 J ? f Q+ . . . 
 (2.3^)aT(c,d) Mn-l)f Q - ^-{(n-l)(Ad 4 )-6c 2 d 2 j^f 0+ i|^/(n-l)(c 6 + d 6 )-15c 2 d 2 (c 2 + d 2 ) ) )^f D 
 
 8h 8 rr / lV 8 .8v O o 2^2, k Ik „ kMJS 
 + 757— U(n-l)(c +d )-28c d (c +d )+70c d jtr 
 
 -2[(n-l)(c 8 + d 8 )- 2 8c 2 a 2 (cVd U )-70oV^] f Q 
 
 + .. . 
 
 (B.35)<C(e) |n(n-l)(n-2)f - ^|-(n-2)(n-7)^f + ^ 2 | r 2 - (n 2 -33n + 122)C7 6 f 
 
 Q O 
 
 ^{(n-2)(n + 13)T3 8 -2(n 2 -129n + 10^^ f ( 
 
 We recall that is the origin, or centre of the square, a (a) includes all points 
 with one coordinate + ah and the rest zero, P(b) has two coordinates each 
 independently + bh and the rest zero, V ( c, d) has one coordinate + ch, another + dh 
 and the rest zero a.nd finally£(e) has three coordinates each independently + eh 
 with the rest zero. 
 
 3. Expansions over a. square are simpler since.!/ f„, ^~ f_ etc. vanish. They can 
 "be obtained by analysis with the detached operators - in particular /vj we proceed 
 to obtain expansions with general terms . 
 
- 3 - 
 
 If F(z) = u + iv is a function of a complex variable z = x + iy then both 
 
 2i 2i 
 
 u and v are harmonic functions satisfying D cj) + D 9 = 0. Likewise, if u is a 
 
 harmonic function, it can he shown that v exists such that u + iv is a function 
 of a complex variable. We then have 
 
 D F = IF' = iD F 
 y n 
 
 D D 
 
 x y 
 
 ■**-< 
 
 ID = -i D 
 x y 
 
 and 
 
 (3-D 
 
 In order to develop expansions we therefore substitute 
 
 (3.2) d x = i'^lr D y = i 1 /^ 
 
 Consider, firstly 
 
 _ /-h s*l . hD -hD hD -hD 
 
 (3.3) J = (2h)* 2 / / f(x,y)dx dy = X (e X -e x )(e y -e y )f 
 
 -h u -h kh D D ^ 
 
 x y 
 
 The operator is 
 
 sinh h®x 3inh hPy s inh i ' ~h^"sinh i ' h % T 
 
 h D D 
 
 x y 
 
 (3.10 
 
 h 2 ^ 2 
 
 _ ,.1/2 .-1/2^, ^ . ,.1/2 .-1/2^^ 
 
 1 cosh (1 / + 1 / )h73 - cosh (1 ' -1 ' )hFT 
 
 = 2 h £ S 2 
 
 1 cosh y2 h -cos >|2 h'fcr 
 
 "2 , 2 _2 
 
 h TJ 
 
 whence 
 
 (1 *\ t I* 2V 4 ^ 2 5 h 8 -JB 2 2r+ V r Lr , . 
 
 ahD -ahD a,bD -ahD 
 
 Likewise/' Z f(x , y ) = (e + e X + e y + e 
 
 = 2 (cosh ah D + cosh ah D ) f_ 
 v x y 
 
 ) tr 
 
 = 2 (cosh 
 
 i" 1 / 2 ahJfr+ cosh l 1 ' 2 ah*5) t r 
 
 (3.6) 
 
 = h cosh 
 
 - Mi - 
 
 ah 
 
 ^<?D ?r 
 
 2"^ 
 
 a h s& m- a, 
 
 $ 
 
 8 U 8 q ^ r v> r ). 
 
 ]f, 
 
- k - 
 
 and 
 
 Z f(x_,y_) = k cosh bh D cosh bh D f_ 
 f p(b) P P X y ° 
 
 = 4 cosh i-^bhJbcosh W 2 bh Sb f 
 
 = 2 (cosh bh N /2~£T+ u bh \/~2^0 f 
 
 , r -, 2 b h <=>.4 2 b h ^.d 
 
 = k [1+ n T3 + 5! W + . 
 
 IT] 
 
 8: 
 
 -2r.kr,kr \ 
 
 b h ^ + ..]f 
 
 V^T 
 
 o 
 
 We shall not use all the expansions given above in the present note, but 
 it seems useful to set out the collected results for future use, 
 
 k. Lattice-point formulae over a, Square . We consider first formulae in two 
 variables, and start with 9 points, putting a = b = 1 and using the sets 0, Qi(l), 
 3(l). We write 
 
 (4.1) J = l/4h = A Q f(0,0) + Z A q f ( Xq , y a ) + Z A p f (Xp,y p ) 
 using (x y ) etc. as typical sets of coordinates. 
 
 Using (3-5) - (3-7), we equal coefficients of Tj f n , r = 0(l)2. 
 
 This gives 
 
 (*.2) 
 
 A A + U + U = 1 
 a 3 
 
 k 
 
 - h A + 16 A = - 
 
 a 3 15 
 + 4 A + 6, A Q = i| 
 
 12 
 
 with correction term C = - ( -k A + 256 A a - -=-=-) — 
 
 OJ p yi -L2 
 
 r^ Sf o 
 
 We obtain the formula 
 
 (fc.3) 
 
 ^ 
 
 -32 
 
 7 
 
 
 -32 
 
 1000 
 
 -32 
 
 
 7 
 
 -32 
 
 7 
 
 ■7 900 
 
 14.1. 4 4-- 19^2 h^ «_L2 - 
 
 with main correction term- ^> ^^T ~ ^ 
 
 This formula is remarkably good. With the example of Part I, we have, 
 
 2 
 writing J ' - h J 
 
 1 /1.2 /d.2 x 
 
 J' = r- / sin x sinh y dx dy = t- (l-cos 1.2)(cosh 1.2 -l) 
 
 Formula, (4.3) gives 
 
 ,12 
 
 = 0.12922 70590 73675 11602 
 
 J T = 0.12922 70590 72834 11029 
 12/ 
 
 withE= -0.-0 84l 00573 and C = +0.0 8^1 01633. 
 
- 5 - 
 5. Five-point formulae . The high precision of (^.3) suggests that formula.e 
 of lessor precision, with fewer points, may be useful. We use the first two 
 of (^.2) and take one of A„, A , A to be zero|(i) A = gives an eight -point 
 formula with relatively poor precision. 
 
 (5.1) 
 
 19 
 
 56 
 
 19 
 
 
 56 
 
 
 
 56 
 
 
 19 
 
 56 
 
 19 
 
 (ii) A 
 
 = 
 
 gives 
 
 
 
 (5.2) 
 
 l 
 
 
 
 l 
 
 
 
 
 
 56 
 
 
 
 
 
 l 
 
 
 
 i 
 
 
 300 
 
 ,8 g 
 with main correction term -k-0 frrTr f 
 
 
 
 6o 
 
 with main correction term - 
 
 32 h 8 _8 
 - 9T ^ f o 
 
 (iii) A. = gives 
 P 
 
 (5.3) 
 
 
 
 -1 
 
 
 
 -1 
 
 19 
 
 -l 
 
 
 
 -l 
 
 
 
 7 15 
 
 with main correction term + — — 
 
 5 9 
 
 r* 8 
 
 
 
 7 8 
 
 We observe tha.t (5.2) and (5«3) combine in the proportions ir= * r-=- to give (^.3).? 
 
 7 k 1J J-5 
 
 though without an error estimated Likewise *■ X (5«2) - -~ >( (5-3) gives (B) of 
 
 lip- 5 1.8 a 
 Note I, and an estimate for the correction, namely - ±±£ £_ *^)° ^ when (x v) is 
 
 harmonic . 
 
 Another combination, that of (5-2) and (5.3) in equal proportions gives a 
 small correction term: 
 
 (5.M 
 
 l 
 
 -k 
 
 1 
 
 
 
 -k 
 
 132 
 
 -k 
 
 
 
 l 
 
 -k 
 
 L 
 
 
 Again k 
 
 K (5.2) 
 
 - 3 
 
 *(5.3 
 
 ) § 
 
 (5.5) 
 
 1 
 
 3 
 
 l 
 
 
 3 
 
 -1 
 
 3 
 
 
 l 
 
 3 
 
 l 
 
 * 120 
 
 with main correction term- p ttt <: W f~ 
 
 5 9» 
 
 r 15 
 
 with main correction term - 
 
 212 
 
 8 
 
 Sr* 8 
 
- 6 - 
 
 Evidently (4.3) is most precise, "but simultaneous use of (5.2) and (5.3) 
 gives an idea, of the precision attained and readily yield the better result if 
 desired. Formula (5«5) might be helpful with desk computing but (5.1) has 
 little to recommend it. 
 
 Numerical results for some of the formulae using the example of £ 4 are as 
 follows : 
 
 Formula. 
 
 Result J' 
 
 
 10 10 / 
 
 E 
 
 
 10 10 K c 
 
 (5.1) 
 
 0.12922 72986 
 
 
 +2395 
 
 
 
 -2396 
 
 (5.2) 
 
 0.12922 70974 
 
 
 + 383 
 
 
 
 - 383 
 
 (5.3) 
 
 0.12922 70255 
 
 
 - 336 
 
 
 
 + 335 
 
 I (B) 
 
 0.12922 71932 
 
 
 +1341 
 
 
 
 +1342 
 
 (5 JO 
 
 0.12922 70615 
 
 
 + 24 
 
 
 
 - 2k 
 
 General n; 
 
 2 
 2n + 1 points. 
 
 We 
 
 consider 
 
 now 
 
 the 
 
 n-dimensio 
 
 n > 3, using lattice points 0, Q:(l), |3(l). In this case the term in J f n is 
 
 o 
 
 relevant, and the 4L f n term will appear in the error, except when n = 3. 
 
 k *J& 
 We equate coefficients of f~, sT f_, u f in the expansions resulting 
 
 from use of (2.1), (2.3l)-(2.33) in (4.1). We obtain 
 
 (6.1) ( A Q + 2n A Q + 2n(n-l) A p = 1 
 
 - k A - 8(n-4) A Q = ~ 
 a v ' 3 15 
 
 + 6 A a + I2(n-l6) Ap = ~ 
 
 o o 
 
 while . -{k A a+ 8(n + 6)A p - §§r*\ + K ^(n-6k)A^ + iff } jjj. s} f Q 
 
 These yield 
 
 ((, o\ A -6ln 2 + 931n + 3780 6ln -496 -6l 
 
 K °' d) A " 3780 a " 3780 A p * T55o 
 
 with „ II98 h 8 ^8 _ 3619 h 8 <±Q _ 
 
 c - -95^ BT* 7 f + -315 BT Xf o 
 
 In particular 
 
 ^ a^ a a 12048 . 626 . 61 
 (6 ' 33) n = 3 A = "75^ A a = " 7560 A P = " 75SO 
 
 (6 34) n - 4 A - i32S§ A St A = - — & 
 
 ^o.j4; n - 4 a q -^q A a ^q A p y^Q 
 
 A a " 7560 
 
 V- 
 
 A 5 ° U 
 
 a " 756o 
 
 V- 
 
- 7 - 
 tti ^^ c « 1382 ° a 382 . 6l 
 
 (6 36) n-6 A =^£ A _._260 A -__4i 
 lD,jDj n _ D A T5S0 A a 75^0 V 75^0 
 
 As a numerical illustration for n = 3 we consider 
 
 1 1 r 1 r 1 C 1 3 5 
 
 J = 7T I = 7T cos jf x cos y cosh f z dx dy dz 
 
 J .1 J-i Jjl 
 
 = ~ sin r sin 1 sinh r = O.98OO827 
 15 4 4 ' 
 
 Formula (6.33) gives J = 0*9799734 WithE= -0.0000109 and C = +0.0000110. 
 
 This result is less spectacular than that of Sk, for these reasons: 
 
 i) In &4, h = 0.6., here h = 1, a.nd the correction term in (4.3) contains a high 
 
 power of h. 
 
 8 12 
 
 ii) The correction term in (6.2) is of order h , that in (4.3) is of order h 
 
 iii) The higher the number of dimensions., the more individual terms there are 
 
 in cD f } JD^ f<> etc. In (4.3) there is only one term in<£) f 3 in (6.33) 
 
 o 
 
 there are 9 i n ^ f* 
 
 iv) The effect of larger interval h is enhanced by the use of the factor 
 
 5 5 
 
 x- y which exceeds unity., in cosh t- Z; this is only partially balanced by 
 
 3 
 
 the factor cos j~ x. 
 
 In spite of these points , the formula (6.2) seems a good one. 
 
 7. Quadrature over a. square^, specially chosen points. Since the expansions 
 of § 3 involve only cross-differences 2) f_, it appears likely that use of 
 sets of diagonal points f3 will be mo."- profitable than attempts to use sets a. 
 It turns out that sets 0, a(a), p(b) and 0,a(a.) both fail to give real values 
 of a if maximum precision is sought. On the other hand,, we can get several 
 formulae making use of any number of sets ^(b^) , -\> = 0(l) r, both with and 
 without the point 0. 
 
 We start first with r sets p ^bp)., without the point 0. We have to find 
 the 2r constants A 3 bj, satisfying the equations 
 r 
 
 (T-1) 1 , k \ b ^ (S " l) = T2 s-l)Os-3) " \-V s " 1(1) 2r 
 
 obtained by substitution of (3.5) — (3. 7) i n 
 (7.2) J = I A ft f (+ b. h, + b. h) 
 
- 8 - 
 
 and equating the coefficients of the first 2r coefficients of - 
 
 powers of k have been cancelled. 
 
 k 
 By familiar arguments, the bi, are rqots of the equation 
 
 2 
 
 Sundry 
 
 (7.3) 
 
 x 
 
 c 
 
 c c - c _ 
 
 r r-1 r-2 
 
 r+1 
 
 '2r-l 
 
 = 
 
 These are the orthogonal polynomials corresponding to the weight function 
 
 1 -3/^- -l/2 
 w(x) = ;r (x ' -x ' ) for range < x < 1. The first two are 
 
 (7.10 
 
 I5x -1=0 
 
 8l9x 2 - k3Qx + 11 = 
 
 (7.5) 
 
 The main correction term is obtained from the next power of S) and yields 
 
 8r x M, 8r <s8r 
 
 C .. (C_ - Z k JL b° r ) 2^ 
 
 2r 
 
 If the point is included, our equations (7»l) are replaced by 
 
 A„ + klA„ = 1 
 
 
 
 3% 
 
 (7.6) 
 
 \>=i' 
 
 £ 4 A„ b 
 
 biL % * = (2s + l)(Wl) = S 
 
 = C , s = 1(1) 2r 
 
 k 
 
 and the b*> are roots of the equation 
 
 2 
 
 1 
 
 X 
 
 X 
 
 c l 
 
 C 2 
 
 C 3 
 
 °2 
 
 C 3 
 
 c k 
 
 °r+l 
 
 °r+2 
 
 °r+3 
 
 'r+1 
 
 'r+2 
 
 '2r 
 
 = 
 
- 9 - 
 
 which are the orthogonal polynomials for the weight function w(x)= ^(x ' -x ' 
 for the range < x < 1. The first two are 
 
 3x - 1 = 
 (7-8) 2 
 
 I7017x - 1365OX + 1745 = 
 
 The main correction term is this turns. 
 
 r o . 2r+l ,4r+2 i, „ 
 
 (T.9) c . (c 2r+1 -£ 4 ^ ^) ^f^ »^ 2 
 
 In each case the coefficients A c may be computed hy standard methods. 
 
 p 
 
 K r 
 
 8. Formulae for r = 1. These ha,ve 4 and 5 points respectively 
 
 8^8 
 
 (8.1) A p = ¥ b = 15 C = 225 "ST"* 
 
 (8.2) A n = 
 
 4 
 
 o 
 
 12 iv 12 
 
 a 1 h q-ia p - 2816 h^ ^ f 
 
 " 5 P " 20 ° " 5 U ~ 122H5 121 J ° 
 
 Written out in full: 
 
 (8.3) j -J 1 -\ (fdr^ir^j^-is-^ir^Hfds-^-is-^M-is^s-^j 
 
 (8.4) j. I 1= | f(o , )^| f (3- 1 A 3 -^),f ( -3-^,3-^) + f(3- lA ,-3-^).f(-3-^-3 1 / i )} 
 
 Asa numerical test use 
 
 t i t 1 r x r 
 
 cos x cosh y dy = sinl sinhl 
 = O.98889 77057 62865 
 
 Formula (8.3) gives 0. 98889 06525 
 
 with E = 0.00000 70533 and C = +0.00000 705+7 
 o^tuL formula (8.4) gives 0. 98889 77062 +1358 
 
 with E,.= +0„0 9 4 78^93 and C = -0. 0^785^3. 
 
 9° Formulae for r = 2 . These have 8 and 9 points respectively 
 
 (9.1) b_ = 0.40316 26030 59346 89754 A Q = 0.22912 30654 28169 97222 
 
 b 2 = 0.84439 75319 23478 74713 Ap = 0.02087 69345 71830 02778 
 
 54592 256 ,16 0.16 
 with mam correction tem—r^nr X TZT n *D f 
 
- 10 - 
 
 (9.2) b Q = A Q = 0.69521 80834 12925 81989 
 
 b n = 0.63205 02078 I8796 99524 A c = 0.06686 42185 46105 38-162 
 
 1 f3]_ 
 
 b 2 = 0.89531 6379I 24106 97730 A^= 0.00993 12606 00663 16340 
 
 ,4.x. +• 16832 v 1024 ,20 q,20 . 
 
 with main correction term w q q a pn , h JO f 
 
 For J = i 
 
 r 1 ri 
 
 cos x cosh y d/ dy 
 
 -1 v/-l 
 
 formula (9.1) gives 0. 98889 77057 62853 38396 
 
 withE= -0.0 13 II71243 and C = +0.0 13 1171555 
 
 while formula (9.2) gives 0. 98889 77057 62865 09647 
 
 with E = +0.0 1 99 and C= 0.0^90 
 
 ■With formula (9«2) we find approximately 
 0.82447 37090 77903 16756 for 
 
 1 r 2 r 2 
 
 -r-?- I cos x cosh y dx dy = sin 2 sinh 2 
 
 = 0.82447 37090 77809 15433 
 with E = 0.0 13 9401323 and C = 0.0 13 94o6250. 
 
 These formulae clearly have high precision, even with considerable values 
 of h. 
 
 10. Quadrature over a. Cube; specially chosen points. The search for such 
 formula.e is more difficult in 3 or more dimensions. It seems that one or more 
 extra available constants are needed in order to obtain real points. We shall 
 not pursue this, but give one simple formula for three dimensions. 
 
 We find nothing convenient by use of points a(a), with or without 0} likewise 
 with (3(b) fails to give real points. We can, however, use 12 points p(b) 
 alone. We have then to satisfy 
 
 2n(n-l)A Q = 1 
 
 4 4 
 8b (4-n)A = — , where n = 3- 
 
 with main correction term. 
 
 l H 21' 6.' ~ " —~'~- " ^ *o 
 
 This yields b = (2/5) ' = 0,79527 07287 67051 A = l/l2 
 
 c = (156 b 6 a + |f) |j. a 6 f - 0.00562 h 6 a 6 f c 
 
- 11 - 
 
 With the example of £6, with integrand cos r- x cos y cosh z 
 (10. 1) gives J = 0.97519 with £ = -0.00489 and C = +0.00^. 
 
 -.6 
 
 The only formula found that allows for the term sj f and ha.s a.n error 
 o 
 
 of order h is (6.33).? which needs 9 points. It is evident that further search 
 is needed. 
 
 This work was supported in part by the Office of 
 Naval Research under Contract Nonr-l83U(27) . 
 

M