Section 5A:

Iterated Gaussian Quadrature over a non-rectangular region

For this section, we modify the 2-dimensional iterated Gaussian rule. Instead of working with a rectangle R = [a,b] x [c,d], we will work with a region S bounded by lines x = a and x = b and by curves y = c(x) and y = d(x).

.

For example, with a = 2, b =5, c(x) = -x and d(x) = x2 we have:

If f(x,y) is continuous on S, then from calculus we know that

I  S f(x, y)dx dy

= [a,b] [c(x),d(x)] f(x, y)dy dx.

We want to get an iterated Gaussian quadrature rule for this kind of region .

For simplicity, we will illustrate the method only for n = 2.

In the x direction we will use Gaussian Quadrature in [a,b]. For example, with n = 2, we get nodes x0, x1, x2 in [a, b] and positive weights w0, w1, w2 :

In the y direction, we apply Gaussian Quadrature for degree 2n+1 on each interval,

[c(x0), d(x0)], [c(x1), d(x1)], [c(x2), d(x2)].

We get points (x0, y00), (x0, y10), (x0, y20) and weights v00, v10, v20 for

interval [c(x0), d(x0)], points (x1, y01), (x1, y11), (x1, y21) and weights v01, v11, v21 for interval [c(x1), d(x1)], and points (x2, y02), (x2, y12), (x2, y22) and weights v02, v12, v22 for interval [c(x2), d(x2)].

We get a grid of (n+1)2 nodes as follows:

The quadrature rule is then:

Q(f) = 0i20j2 wivij f(xi, yij)

For a general n (replacing 2 by n), we get

Q(f) = 0in0jn wivij f(xi, yij),

This method can be computed in Mathematica as follows:

The code is a revised version of the function G2dim from Section 4A. This function uses the method Gaussian (Section 3A), which finds Gaussian points

x[0], …, x[n] and weights weights[1],….,weights[n+1]. Gaussian is used once in the x-direction, for the interval [a,b], and n+1 times in the y-direction, on each y-interval

[c(x[i]), d(x[i])] (0 i  n).

The final answer is stored in variable II. In the next section, we illustrate this method

of approximating integrals.

There is one important difference between this method and the method in Section 4A for rectangular regions. The method in Section 4A is exact if the degrees of f(x,y) in x and y separately are both  2n+1. This is not always true for the new method. The problem is that even if the degree of f(x,y) in x is at most 2n+1, the degree of

G(x)  [c(x),d(x)] f(x, y) dy

in x may be greater than 2n+1. Now,

I S f(x, y)dx dy

= [a,b] [c(x),d(x)] f(x, y)dy dx

= [a,b] G(x) dx.

When we use Gaussian quadrature in the x-direction to estimate the last integral,

we cannot expect exact results because the degree of G(x) may be greater than 2n+1

(see Example 3 below).

Section 5B:

Experiments with 2-dimensional iterated Gaussian Quadrature over a non-rectangular region

The region we consider is bounded by by the lines x = a and x = b, and by the curves y = c(x) and y = d(x). The method is described in Section 5A. In these experiments we are trying to find the smallest n which gives a Gaussian relative error less than 1%

In the first example, we use c(x) = c and d(x) = d to show that the method agrees with the previous method for rectangular regions (Section 4A-4B). Throughout these experiments, we use 100-digit precision. Here is the first example.

Example 1:

f(x,y) = x7y7 over the regions 1 x  2 , 3  y  5 (c(x) = 3, d(x) = 5)

We see that n = 2 satisfies our objective. Since the degrees in x and y separately are 7

( = (2 x 3) +1), we expect exact results for n = 3.

Mathematica / Iterated Gaussian / Gaussian
Relative Error
N = 1 / 1.53026*106 / 1.48531*106 / 0.0293733
N = 2 / 1.53003*106 / 0.000144306
N = 3 / 1.53026*106 / 2.74423 x10-10

Due to roundoff error, we did not get exact results for n = 3, but we did get a relative error very close to zero.

Example 2:

f(x,y) = 1 over the region –1  x  1, -Sqrt(1-x2)  y  Sqrt(1-x2)

This estimates the volume of a right circular cylinder with radius one and height one, so the exact answer is .

We see that n = 3 satisfies our objective.

Mathematica / Iterated Gaussian / Gaussian
Relative Error
N = 1 / 3.14159 / 3.26599 / 0.0395957
N = 2 / 3.18323 / 0.013255
N = 3 / 3.16056 / 0.00603593
N = 6 / 3.14556 / 0.0012641
N = 8 / 3.14353 / 0.000618037
N = 20 / 3.14176 / 0.000052877
N = 22 / 3.14172 / 0.0000404815

We tried to find the smallest n so that the Gaussian approximation agreed with  to 6 places, but Mathematica would not execute the code for any n  23. The problem was that Mathematica could not correctly invert the matrices in the Gaussian code, because for large values of n, these matrices are very badly conditioned.

Example 3:

f(x,y) = xy2 over the region 2  x  5, -x  y  x2 (see picture in Section 5A)

n = 2 satisfied our objective of an error under 1%.

Mathematica / Iterated Gaussian / Gaussian
Relative Error
N = 1 / 16471.6 / 15766.3 / 0.0428169
N = 2
N = 3 / 16465.2
16471.6 / 0.000387258
0

Notice that in this example the degrees of f in x and y separately are both less than

3 (= 2  1 + 1), but we do not get an exact result for n = 1. This is an example of the problem that is mentioned at the end of Section 5A. In this case,

G(x)  [c(x),d(x)] f(x, y) dy

= xy3/3 |[-x,x2] = (x4 + x7)/3, so deg G(x) = 7. So we cannot expect Gaussian quadrature

with n =1 to give exact results for the integral of G. Since 7 = 2  3 + 1, we do expect

exact results for n = 3, and this is what the table shows.

Example 4:

f(x,y) = e(x-y) over the regions 2  x  5, -x  y  x2

Up to n = 5, the Gaussian relative error was greater than 1%. However, n = 6 gave a small error.

Mathematica / Iterated Gaussian / Gaussian
Relative Error
N = 1 / 10985.9 / 955.728 / 0.913004
N = 2 / 10188.6 / 0.0725768
N = 3 / 10188.6 / 0.0725768
N = 4 / 10188.6 / 0.0725768
N = 5 / 10798.8 / 0.017034
N = 6 / 10950.8 / 0.003196
N = 7 / 10980.5 / 0.000490384

Example 5:

f(x,y) = Sqrt[1-(x2 + y2)] over the region –1  x  1 , -Sqrt[1-x2]  y  Sqrt[1-x2]

Mathematica / Iterated Gaussian / Gaussian
Relative Error
N = 1 / 2.0944 / 2.17732 / 0.0395957
N = 2 / 2.12216 / 0.013255
N = 3 / 2.10704 / 0.00603592
N = 4 / 2.10121 / 0.00325313

We are estimating the volume of the upper unit hemisphere.

The volume of a hemisphere is V = 4/3r3. Since we are only dealing with the upper half, the volume is V = 2/3r3. In example 5, r = 1. So the volume of the hemisphere is V = 2/3 which is approximately 2.0944. We get a relative error under 1% for n = 3.