Subsections

• B.3.1 Gaussian Quadrature
• B.3.2 Integral Transformation
• B.3.3 Legendre Polynomials
• B.3.4 Numerical Quadrature Library

B.3 Volume Integral

The numerical evaluation of a volume integral over an arbitrary tetrahedron should, if possible, be widely avoided. Under different circumstances this rule of thumb cannot be fulfilled and the software engineer rely on a computational evaluation. In this case it is highly recommended to use existing code or libraries as described in Section B.3.4. Nevertheless, the following gives a rough sketch of a possible numerical volume integral evaluation.

1 Theorem ($3$ Dimensional Change of Variables Theorem)   Let $V$ and $W$ be a bound open sets in $\mathbb{R}^3$ and $\mathbf{R} : V \rightarrow W (\mathbf{R}(u,v,w) = (X(u,v,w), Y(u,v,w), Z(u,v,w))$ be a $C^1$ map which is one to one and onto. Suppose that $f: W \rightarrow \mathbb{R}$ is a continuous function, then

$$\iiint\limits_W f(x,y,z) dx dy dz = \iiint\limits_V f(\mathbf{R}(u,v,w))J(u,v,w) du dv dw,$$ (B.11)

where

$$J(u,v,w) = \Bigg \vert \det \begin{bmatrix}\mathbf{R}_u (u,v,w) \\ \mathbf{R}_v (u,v,w) \\ \mathbf{R}_w (u,v,w) \end{bmatrix}\Bigg \vert \notag$$ (B.12)

$$ = \Bigg \vert \det \begin{bmatrix}X_u (u,v,w) & Y_u (u,v,w) & Z_u (... ...v (u,v,w) \\ X_w (u,v,w) & Y_w (u,v,w) & Z_w (u,v,w) \end{bmatrix} \Bigg \vert$$ (B.13)

is called the Jacobian of the transformation $\mathbf{R}$ .

Remark 6   The content of this theorem may be summarized by the equation

$$dx dy dz = \vert \det \mathbf{R}' (u,v,w) \vert \; du dv dw$$ (B.14)

where

$$\mathbf{R}' (u,v,w) \equiv \begin{bmatrix}\mathbf{R}_u (u,v,w) \ \mathbf{R}_v (u,v,w) \ \mathbf{R}_w (u,v,w) \end{bmatrix}.$$ (B.15)

The volume integral over an arbitrary tetrahedron mapped to a standard $3$ -simplex can now be written as

$$\int_\mathcal{V} f(x,y,z)\; dV = \int_{0}^{1} \int_{0}^{1-\xi } \int_{0}^{1-\xi - \eta }f(\mathbf{... ...vert \det \mathbf{R}' (\xi ,\eta ,\zeta ) \vert \; d\xi d\eta d\zeta \notag$$ (B.16)

$$ = \int_{0}^{1} \int_{0}^{1-\xi } \int_{0}^{1-\xi - \eta }f(\mathbf{R}(\xi ,\eta ,\zeta )) \; \vert \det \mathbf{J} \vert \; d\xi d\eta d\zeta$$ (B.17)

where $\mathbf{J}$ is the Jacobian matrix given in (3.5).

Numerical integration methods can generally be described as combining evaluations of the integrand to get an approximation of the integral. An important part of the analysis of any numerical integration method is to study the behavior of the approximation error as a function of the number of integrand evaluations. A method which yields a small error for a small number of evaluations is usually considered superior. Reducing the number of evaluations of the integrand reduces the number of arithmetic operations involved, and therefore reduces the total round-off error. Also, each evaluation consumes time and the integrand may be arbitrarily complicated.

B.3.1 Gaussian Quadrature

A large class of quadrature rules can be derived by constructing interpolating functions which are easy to integrate. Typically these interpolating functions are polynomials. The Newton-Cotes formulas (named after Isaac Newton and Roger Cotes) are a group of formulas for numerical integration. It is assumed that the value of a function $f$ is known at equally spaced points $x_i$ , for $i=0,\ldots,n$ . If the evaluation points are not assumed to be equally spaced, another class of formulas, called Gausssian quadrature, can be derived. With the Gaussian quadrature rule an exact result for polynomials of degree $2n+1$ , by a suitable choice of the $n$ points $x_i$ and weights $w_i$ can be achieved. The domain of integration for such a rule is conventionally taken as $[a,b]=[-1,1]$ , so the rule is stated as

$$\int_{-1}^1 f(x) dx \approx \sum_{i=1}^n w_i f(x_i).$$ (B.18)

It can be shown that the evaluation points are just the roots of a polynomial belonging to a class of orthogonal polynomials [136].

Table B.1: Some polynomials for general Gaussian quadrature.

Interval	$\omega(x)$	Orthogonal polynomials
$[-1,1]$	$1$	Legendre polynomials
$[-1,1]$	$\frac{1}{\sqrt{1-x^2}}$	Chebyshev polynomials
$[0,\infty)$	$e^{-x}$	Laguerre polynomials
$(-\infty,\infty)$	$e^{-x^2}$	Hermite polynomials

Table B.1 gives an overview of some polynomials used for the more general Gaussian quadrature where by introducing a weight function $\omega(x)$ into the integrand and allowing interval others than $[-1,1]$ the problem is given by

$$\int_{a}^b \omega (x) f(x) dx.$$ (B.19)

If an integral within the interval $[a,b]$ must be changed into an integral with the limits $[-1,1]$ , the following transformation can be applied:

$$\int_a^b f(t) dt = \frac{b-a}{2} \int_{-1}^1 f\left(\frac{b-a}{2}x + \frac{a+b}{2}\right) dx$$ (B.20)

obtaining the following approximation

$$\frac{b-a}{2} \sum_{i=1}^n w_i f\left(\frac{b-a}{2}x_i + \frac{a+b}{2}\right).$$ (B.21)

The error of a Gaussian quadrature rule can be stated as follows. For an integrand with $2n$ continuous derivatives,

$$\int_a^b \omega (x) f(x) dx - \sum_{i=1}^n w_i f(x_i) = \frac{f^{(2n)}(\xi)}{(2n)!} \Vert p_n\Vert^2$$ (B.22)

for some $\xi$ in $(a,b)$ , where $p_n$ is the orthogonal polynomial of order $n$  [137].

B.3.2 Integral Transformation

For a quadrature approximation of the volume integral given in formula (B.15), one has to take into account that the points of evaluation reflect the nature of the integration volume. For example, for an integration along the $\zeta$ coordinate, the according differential volume slice of the standard $3$ -simplex gets smaller and smaller with increasing $\zeta$ values. A common way of integration is to reduce the problem from a simplex to an integration over a cube [138,139]. A general extension of this approach regarding the standard $n$ -simplex is given in [140]. Based on the work given in [141] the polynomial class of Legendre polynomials is chosen for the quadrature of a standard $3$ -simplex mapped on a cube.

An integral over the standard $3$ -simplex may be written as

$$\int_{\mathcal{V}_{S3}} f(x,y,z)\; dV = \int_{0}^{1} \int_{0}^{1-x} \int_{0}^{1- x -y} f(x,y,z)\; dz dy dx.$$ (B.23)

With the following transformation,

$$x = x$$

$$y = \hat{y} (1-x)$$ (B.24)

$$z = \hat{z}(1-x-y) = \hat{z}(1-\hat{y})(1-x)$$

the limits are transformed to 0 and $1$ . Since the Jacobian of the above transformation is given by

$$J = (1-x)^2 (1-\hat{y}),$$ (B.25)

the volume integral can be rewritten as

$$\int_{0}^{1} \int_{0}^{1-x} \int_{0}^{1- x -y} f(x,y,z)\; dz dy... ...0}^{1} (1-x)^2(1-\hat{y}) \hat{f}(x,\hat{y},\hat{z})\; d\hat{z} d\hat{y} dx$$ (B.26)

and finally a Gaussian quadrature based on Legendre polynomials can now be used to solve formula (B.24).

B.3.3 Legendre Polynomials

Legendre polynomials are solutions of Legendre's differential equation defined as

$$\frac{d}{dx} \left[ (1-x^2) \frac{d}{dx} P(x) \right] + n(n+1)P(x) = 0.$$ (B.27)

The Legendre differential equation can be solved using the standard power series method. The solution is finite (i.e. the series converge) provided for $\vert x\vert < 1$ . Furthermore, it is finite at $x = \pm 1$ provided $n$ is a non-negative integer, $n=0,1,2,\ldots$ . Than the solutions form an orthogonal polynomial sequence called the Legendre polynomials.

Each Legendre polynomial $P_n(x)$ has degree $n$ and is expressed via

$$P_n(x) = (2^n n!)^{-1} \frac{d^n}{dx^n } \left[ (x^2 -1)^n \right].$$ (B.28)

Their orthogonality is defined by the expression

$$\int_{-1}^{1} P_m(x) P_n(x) dx = \begin{cases}0& \text{if $m \neq n$} \ \frac{2}{2n+1}& \text{if $m = n$}. \end{cases}$$ (B.29)

The first few Legendre polynomials are

$$P_0(x) = 1$$

$$P_1(x) = x$$

$$P_2(x) = \frac{1}{2}(3x^2-1)$$ (B.30)

$$P_3(x) = \frac{1}{2}(5x^3-3x)$$

$$P_4(x) = \frac{1}{8}(35x^4-30x^2+3)$$

and their roots $(r_0 \ldots r_n)$ can be found numerically by well-known methods [142]. According to formula (B.16), for the weights $(w_0 \ldots w_n)$ one has to solve the linear system given as follows:

$$\left( \begin{array}{cccc} 1 & 1 & \cdots & 1 \ r_0^1 & r_1^1 & ... ...}^{1} \xi^2 d\xi \ \vdots \ \int_{-1}^{1} \xi^n d\xi \end{array} \right),$$ (B.31)

which allows an exact integration of polynomials up to degree $2(n-1)+1$ , where $n$ is the degree of Legendre polynomials, i.e. with $P_3(x)$ a polynomial of degree $5$ can be integrated exactly. Table B.2 gives a numerical solution of equation (B.29) for the first few Legendre polynomials up to $P_4(x)$ .

Table B.2: Roots at the Legendre polynomials up to the degree of $4$ . The weights give the values for the Gaussian quadrature, with respect to formula (B.16).

	roots	weights
	$-0.577\;350\ldots$	$1.000\;000\ldots$
$P_2(x)$	$0.577\;350\ldots$	$1.000\;000\ldots$
	$0.774\;596\ldots$	$0.555\;555\ldots$
$P_3(x)$	$0.000\;000\ldots$	$0.888\;888\ldots$
	$-0.774\;596\ldots$	$0.555\;555\ldots$
	$-0.861\;136\ldots$	$0.347\;855\ldots$
	$-0.339\;981\ldots$	$0.652\;145\ldots$
$P_4(x)$	$0.339\;981\ldots$	$0.652\;145\ldots$
	$0.861\;136\ldots$	$0.347\;855\ldots$

B.3.4 Numerical Quadrature Library

Due to the wide spectrum of numerical integration problems, approaches are multifaceted. Several libraries (under different licenses) are available offering more or less out of the box solutions for the problems depicted in this chapter. To give a short indication for further investigations, two libraries are presented.

For a direct evaluation of integrals over the standard $n$ -simplex the Numerical Algorithm Group (NAG) [143] library routine D01PAF can be used. This routine returns a sequence of approximations to the integral of a function over a multi-dimensional simplex, together with an error estimate for the last approximation. This method is primary based on the work given in [144,145]. The computational costs depend mostly on the number of function evaluations, which in turn depends mostly on the highest computed order of the integral approximation $MaxOrd$ and the number of dimensions $Ndim$ for the integral. The number of function calls can be estimated by

$$\frac{(MaxOrd + Ndim)!}{(MaxOrd-1)!(Ndim+1)!}.$$ (B.32)

which gives a rough idea about the complexity of the algorithm.

Another very popular library which offers several different numerical quadrature algorithms, mostly based on polynomial integration schemes, is the so-called GNU Scientific Library (GSL) [146]. This library is a collection of routines for numerical computing. The routines have been written in C and present a modern Applications Programming Interface (API) for C programmers, allowing wrappers to be written for high level languages. The source code is distributed under the GNU General Public License.