Approximation of Definite Integrals Over the Volume of the Ball

TL;DR

This paper introduces a novel RBF-FD inspired method for efficiently computing definite integrals over 3D ball volumes using arbitrarily scattered nodes, avoiding uniformity constraints and achieving high accuracy with optimal computational complexity.

Contribution

The paper presents a new RBF-FD based algorithm that computes quadrature weights for scattered nodes in three dimensions without uniformity restrictions, with improved efficiency.

Findings

Computes quadrature weights in O(N log N) operations.

Achieves high-order accuracy for scattered nodes.

Removes the need for uniform node distribution.

Abstract

A Radial Basis Function Generated Finite-Differences (RBF-FD) inspired technique for evaluating definite integrals over the volume of the ball in three dimensions is described. Such methods are necessary in many areas of Applied Mathematics, Mathematical Physics and many other application areas. Previous approaches needed restrictive uniformity in the node set, which the algorithm presented here does not require. By using RBF-FD approach, the proposed algorithm computes quadrature weights for $N$ arbitrarily scattered nodes in only $O(N\mbox{ log}N)$ operations with high orders of accuracy.

Equations

$$\underset{\Omega }{\iiint }f(\mathbf{x})dV,$$

$$\underset{\Omega }{\iiint }f(\mathbf{x})dV=\sum _{k\in {\mathcal{K}}_I}\underset{t_k}{\iiint }f(\mathbf{x})dV+\sum _{k\in {\mathcal{K}}_S}\left(\underset{t_k}{\iiint }f(\mathbf{x})dV+\underset{s_k}{\iiint }f(\mathbf{x})dV\right).$$

$$\underset{t_k}{\iiint }f(\mathbf{x})dV$$

$$\varphi \left({\Vert \mathbf{x}-{\mathbf{x}}_{k,j}\Vert }_2\right),j=1,2,\dots ,n$$

$$s(\mathbf{x}):=\sum _{j=1}^n{c}_{k,j}^{\text{mbox}RBF}\varphi \left({\Vert \mathbf{x}-{\mathbf{x}}_{k,j}\Vert }_2\right)+\sum _{l=1}^M{c}_{k,l}^p{\pi }_l(\mathbf{x}),$$

$$\underset{t_k}{\iiint }f(\mathbf{x})dV\approx \sum _{j=1}^n{w}_{k,j}f({\mathbf{x}}_{k,j})$$

$$\underset{t_k}{\iiint }f(\mathbf{x})dV+\underset{s_k}{\iiint }f(\mathbf{x})dV\approx \sum _{j=1}^n{w}_{k,j}f({\mathbf{x}}_{k,j})$$

$$\displaystyle A_{k}=\left[\begin{array}[]{cc}\Phi_{k}^{T}&P_{k}\\ P_{k}^{T}&0\end{array}\right].$$

$${\Phi }_{k,ij}=\varphi \left(\Vert {\mathbf{x}}_{k,i}-{\mathbf{x}}_{k,j}\Vert \right),\text{mbox}fori,j=1,2,\dots ,n.$$

$$P_{k,il}={\pi }_l({\mathbf{x}}_{k,i}),\text{mbox}fori=1,2,\dots ,n\text{mbox}andl=1,2,\dots ,M.$$

$$\displaystyle I_{k,j}=\left\{\begin{array}[]{cc}\iiint\limits_{t_{k}}\phi\left(\left\lVert\mathbf{x}-\mathbf{x}_{k,j}\right\rVert\right)dV&j=1,2,\ldots,n\\ \iiint\limits_{t_{k}}\pi_{j-n}(\mathbf{x})dV&j=n+1,n+2,\ldots,n+M\end{array}\right..$$

$$\displaystyle I_{k,j}=\left\{\begin{array}[]{cc}\iiint\limits_{t_{k}}\phi\left(\left\lVert\mathbf{x}-\mathbf{x}_{k,j}\right\rVert\right)dV+\iiint\limits_{s_{k}}\phi\left(\left\lVert\mathbf{x}-\mathbf{x}_{k,j}\right\rVert\right)dV&j=1,2,\ldots,n\\ \iiint\limits_{t_{k}}\pi_{j-n}(\mathbf{x})dV+\iiint\limits_{s_{k}}\pi_{j-n}(\mathbf{x})dV&j=n+1,\ldots,n+M\end{array}\right..$$

$${\mathbf{n}}_{{\mathbf{a}}_k{\mathbf{b}}_k{\mathbf{c}}_k}:=\frac{({\mathbf{b}}_k-{\mathbf{a}}_k)\times ({\mathbf{c}}_k-{\mathbf{a}}_k)}{{\Vert ({\mathbf{b}}_k-{\mathbf{a}}_k)\times ({\mathbf{c}}_k-{\mathbf{a}}_k)\Vert }_2},$$

$${\mathbf{e}}_{k,j}={\mathbf{x}}_{k,j}+\left[({\mathbf{a}}_k-{\mathbf{x}}_{k,j})\cdot {\mathbf{n}}_{{\mathbf{a}}_k{\mathbf{b}}_k{\mathbf{c}}_k}\right]{\mathbf{n}}_{{\mathbf{a}}_k{\mathbf{b}}_k{\mathbf{c}}_k},$$

$$\underset{t_k}{\iiint }\varphi \left({\Vert \mathbf{x}-{\mathbf{x}}_{k,j}\Vert }_2\right)dV$$

$$\displaystyle\Bigg{\{}\mbox{sign}\left(\left(\mathbf{x}_{k,j}-\mathbf{e}_{k,j}\right)\cdot\mathbf{n}_{\mathbf{a}_{k}\mathbf{b}_{k}\mathbf{c}_{k}}\right)\iiint\limits_{t_{\mathbf{x}_{k,j}\mathbf{a}_{k}\mathbf{b}_{k}\mathbf{c}_{k}}}\phi\left(\left\lVert\mathbf{x}-\mathbf{x}_{k,j}\right\rVert_{2}\right)+\cdots$$

$$\text{mbox}sign\left(\left({\mathbf{x}}_{k,j}-{\mathbf{f}}_{k,j}\right)\cdot {\mathbf{n}}_{{\mathbf{a}}_k{\mathbf{d}}_k{\mathbf{b}}_k}\right)\underset{t_{{\mathbf{x}}_{k,j}{\mathbf{a}}_k{\mathbf{d}}_k{\mathbf{b}}_k}}{\iiint }\varphi \left({\Vert \mathbf{x}-{\mathbf{x}}_{k,j}\Vert }_2\right)+\cdots$$

$$\text{mbox}sign\left(\left({\mathbf{x}}_{k,j}-{\mathbf{g}}_{k,j}\right)\cdot {\mathbf{n}}_{{\mathbf{a}}_k{\mathbf{c}}_k{\mathbf{d}}_k}\right)\underset{t_{{\mathbf{x}}_{k,j}{\mathbf{a}}_k{\mathbf{c}}_k{\mathbf{d}}_k}}{\iiint }\varphi \left({\Vert \mathbf{x}-{\mathbf{x}}_{k,j}\Vert }_2\right)+\cdots$$

$$\displaystyle\mbox{sign}\left(\left(\mathbf{x}_{k,j}-\mathbf{h}_{k,j}\right)\cdot\mathbf{n}_{\mathbf{b}_{k}\mathbf{d}_{k}\mathbf{c}_{k}}\right)\iiint\limits_{t_{\mathbf{x}_{k,j}\mathbf{b}_{k}\mathbf{d}_{k}\mathbf{c}_{k}}}\phi\left(\left\lVert\mathbf{x}-\mathbf{x}_{k,j}\right\rVert_{2}\right)\Bigg{\}}.$$

$$\underset{t_{{\mathbf{x}}_{k,j}{\mathbf{a}}_k{\mathbf{b}}_k{\mathbf{c}}_k}}{\iiint }\varphi \left({\Vert \mathbf{x}-{\mathbf{x}}_{k,j}\Vert }_2\right)dV$$

$$\mathbf{x}(\sigma ,{\lambda }_1,{\lambda }_2)={\mathbf{x}}_{k,j}+\sigma ({\lambda }_1{\mathbf{a}}_k+{\lambda }_2{\mathbf{b}}_k+(1-{\lambda }_1-{\lambda }_2){\mathbf{c}}_k-{\mathbf{x}}_{k,j}).$$

$$\underset{0}{\overset{1}{\int }}$$

$${\sigma }_{{\mathbf{a}}_k{\mathbf{b}}_k{\mathbf{c}}_k}=\frac{\left(\frac{1}{3}({\mathbf{a}}_k+{\mathbf{b}}_k+{\mathbf{c}}_k)-{\mathbf{x}}_{k,j}\right)\cdot {\mathbf{n}}_{{\mathbf{a}}_k{\mathbf{b}}_k{\mathbf{c}}_k}}{\left({\mathbf{c}}_k-{\mathbf{x}}_{k,j}\right)\cdot {\mathbf{n}}_{{\mathbf{a}}_k{\mathbf{b}}_k{\mathbf{c}}_k}}$$

$${\mathbf{x}}^{\prime }(\lambda ,\mu )=(1-\lambda ){\mathbf{a}}_k+\lambda \left((1-\mu ){\mathbf{b}}_k+\mu {\mathbf{c}}_k\right),\text{mbox}0\le \lambda \le 1\text{mbox}and0\le \mu \le 1,$$

$$\mathbf{x}(\lambda ,\mu ,\sigma )=\left(1+\frac{\sigma }{{\Vert {\mathbf{x}}^{\prime }(\lambda ,\mu )\Vert }_2}\right){\mathbf{x}}^{\prime }(\lambda ,\mu ),$$

$$\underset{s_k}{\iiint }\varphi \left(\Vert \mathbf{x}-{\mathbf{x}}_{k,j}\Vert \right)dV=\underset{0}{\overset{1}{\int }}\underset{0}{\overset{1}{\int }}\underset{0}{\overset{\rho -{\Vert {\mathbf{x}}^{\prime }(\lambda ,\mu )\Vert }_2}{\int }}\varphi \left(\Vert \mathbf{x}(\lambda ,\mu ,\sigma )-{\mathbf{x}}_{k,j}\Vert \right)\left|J(\sigma ,\lambda ,\mu )\right|d\sigma d\lambda d\mu ,$$

$$\underset{\Omega }{\iiint }f(\mathbf{x})dV\approx \sum _{k=1}^K\sum _{j=1}^n{w}_{k,j}f({\mathbf{x}}_{k,j}).$$

$$\underset{\Omega }{\iiint }f(\mathbf{x})dV\approx \sum _{i=1}^N{W}_{i}f({\mathbf{x}}_i).$$

$$f_1(x,y,z)=\sum _{\alpha =0}^{30}\sum _{\beta =0}^{\alpha }\sum _{\gamma =0}^{\alpha -\beta }{a}_{\alpha \beta \gamma }{x}^{\alpha -\beta -\gamma }{y}^{\beta }{z}^{\gamma }.$$

$$\sum _{\alpha =0}^{30}\sum _{\beta =0}^{\alpha }\sum _{\gamma =0}^{\alpha -\beta }{a}_{\alpha \beta \gamma }\frac{{\rho }^{\alpha +3}}{8\Gamma \left(\frac{\alpha +5}{2}\right)}$$

$$\Gamma \left(\frac{\beta +1}{2}\right)\Gamma \left(\frac{\gamma +1}{2}\right)\Gamma \left(\frac{\alpha -\beta -\gamma +1}{2}\right)]$$

$$f_2(x,y,z)=\exp \left(-10\left((x-x_s{)}^2+(y-y_s{)}^2+(z-z_s{)}^2\right)\right)$$

$$(x_s,y_s,z_s)=(0.047056440432708,0.071766893999009,0.118950756342700)$$

Full text

{AMS}

68Q25, 65R99

1 Introduction

This article is concerned with the development of a method for the approximate evaluation of definite integrals over the volume of the ball in $\mathbb{R}^{3}$ of radius $\rho>0$. That is, consider volume integration over the domain $\Omega=\left\{\mathbf{x}\in\mathbb{R}^{3}|\left\lVert\mathbf{x}-\mathbf{x}_{0}\right\rVert_{2}\leq\rho\right\}$, with $\mathbf{x}_{0}\in\mathbb{R}^{3}$. Applications of this common problem in mathematics abound from the more specific scenarios of estimating the volumes of wells containing hydrocarbons beneath the Earth’s surface [19] and recovering important diagnostic information in, e.g., Thermoacoustic or Photoacoustic Tomogrpahy [9] to some more generic problems in potential theory, magnetism and other concepts in mathematical physics [11].

The approximation of the values of definite integrals (quadrature when considering integration over an interval, or quadrature/cubature when considering integration over domains in two or more dimensions) is a rich topic dating back many centuries to early attempts to measure, for instance, the area of the circle [7]. Since that time much research has been devoted to developing sophisticated and accurate techniques for estimating the values of integrals over intervals, areas, and volumes. There are many texts devoted to summarizing such methods, see, for example [20, 8, 10, 7].

In the simplest case, the rules for quadrature over intervals are often constructed by replacing the integrand with a polynomial interpolant or a polynomial approximation and then integrating, or by enforcing that a rule be exact on a particular class of functions (like polynomials up to a particular degree). This is how to arrive at, for instance, Newton-Cotes or Gaussian Quadrature rules, respectively. These rules for integration over intervals are then often leveraged in the context of iterated integrals by employing one-dimensional quadrature rules for each variable in turn, leading to so-called Cartesian product formulas. Such formulas require structure in the node sets–spacing between nodes that is uniform or tied to the roots of orthogonal polynomials–across each variable of integration. This requirement may be impractical or require an additional interpolation between the structured node set and locations where the integrand is specified.

To overcome these requirements on the structure of the node set, the value of the integral can be approximated utilizing concepts from number theory or through pseudorandom sampling of the integrand as in Monte-Carlo techniques [8]. More common, however, is the alternative method of constructing quadrature rules over areas and volumes by replacing the integrand with a now multivariate polynomial interpolant or approximation and then integrating. It is often the case that the basis set used for interpolation does not depend on the locations of the nodes (e.g. multivariate polynomials). Such basis sets suffer from a question of the existence and uniqueness of an interpolant on unstructured node sets [13, 12], which has prompted the use of Radial Basis Functions (RBFs) in the approximation of the integrand.

The proposed quadrature technique was developed to compute weights at whatever node locations are specified by the user. This is because in applications, numerical quadrature is usually a follow-up to some other task (such as collecting data, or numerically solving PDEs), making it impractical to require node locations that are specific to the quadrature method. The present algorithm is therefore designed to find the quadrature weights given a node set defined by the application. Further, the proposed algorithm allows for node sets featuring spatially varying separation when increased node density is needed to capture fine structure in the integrand.

The numerical method described in this paper is a generalization of RBF-FD (radial basis function-generated finite differences). This approach has so far mostly been used to approximate partial derivatives, with the key difference to regular finite differences that the node points no longer need to be grid-based (in particular, Cartesian node layouts are now known to be less-than-optimal [4]. For surveys of RBFs and of RBF-FD methods as these are applied to PDEs, see [6, 5]. Further, RBFs have been used successfully to construct quadrature rules for an interval in one-dimension, bounded domains in two-dimensions and, more specifically, integrals over bounded two-dimensional (piecewise-)smooth surfaces embedded in three-dimensions [16, 18, 17].

The following Section 2 describes the present quadrature method. Section 3 describes some test examples, with illustrations of convergence rates and computational costs. Finally, section 4 outlines some conclusions. A Matlab implementation of the method is available at Matlab Central’s File Exchange [15].

2 Description of the key steps in the algorithm

Consider evaluating

[TABLE]

where $\Omega=\left\{\mathbf{x}\in\mathbb{R}^{3}:\left\lVert\mathbf{x}-\mathbf{x}_{0}\right\rVert_{2}\leq\rho\right\}$.

Similar to the work presented in [16, 18, 17] for surface integrals, the proposed algorithm can be described in four steps:

Decompose the domain of integration into $K\in\mathbb{Z}^{+}$ subdomains. 2. 2.

On the $k^{\mbox{th}}$ subdomain construct an interpolant of the integrand. 3. 3.

Integrate the interpolant of the integrand to determine weights for integrating a function $f$ over the $k^{\mbox{th}}$ subdomain. 4. 4.

Combine the weights for the integrals over the $K$ subdomains to obtain a weight set for approximating the volume integral of $f$ over $\Omega$.

Each of these steps is described in greater detail in what follows.

2.1 Step 1: Decompose the domain of integration

Consider first a set $\mathcal{S}_{N}=\left\{\mathbf{x}_{i}\right\}_{i=1}^{N}$ of $N$ points in $\Omega$, with a subset exactly on the boundary surface. On this set of points construct a tessellation $T=\left\{t_{k}\right\}_{k=1}^{K}$ (via Delaunay tessellation or some other algorithm) of $K$ tetrahedra. These tetrahedra encompass the bulk of the volume of $\Omega$ but not the entire volume.

Let $\mathcal{K}_{S}\subset\{1,2,\ldots,K\}$ be the set of indices such that if $k\in\mathcal{K}_{S}$ then the tetrahedron $t_{k}$ has a face that is not shared with any of the other tetrahedra. This face has all three vertices on the surface of $\Omega$ and, unless the surface is planar, there is a sliver of volume, $s_{k}$, between the unshared face and the spherical surface bounding the ball that must be accounted for in decomposing the volume. Conversely, let $\mathcal{K}_{I}=\left\{1,2,\ldots,K\right\}\backslash\mathcal{K}_{S}$ be the indices of tetrahedra that do not have a face with three vertices on the surface. With these definitions (1) can be decomposed as

[TABLE]

2.2 Step 2: Construct an interpolant of the integrand

For each tetrahedron in $T$ define the sets $\mathcal{N}_{k}=\left\{\mathbf{x}_{k,j}\right\}_{j=1}^{n}$ to be the $n$ points in $\mathcal{S}_{N}$ nearest to the midpoint of $t_{k}$ (the average of the vertices of $t_{k}$). Then for an individual $t_{k}$ the integral

[TABLE]

is evaluated by first approximating $f(\mathbf{x})$ by an RBF interpolant, with interpolation points from the set $\mathcal{N}_{k}$, and then integrating the interpolant. Often the RBF interpolant is a linear combination of (conditionally-) positive definite RBFs,

[TABLE]

and augmented by multivariate polynomial terms. If $k\in\mathcal{K}_{S}$, then the same interpolant and set of nodes is used for approximating the integrand over $s_{k}$. Define $\{\pi_{l}(\mathbf{x})\}_{l=1}^{M}$, with $M=\frac{(m+1)(m+2)(m+3)}{6}$, to be the set of all of the trivariate polynomial terms up to degree $m$. The interpolant is constructed as

[TABLE]

where $c_{k,1}^{RBF},\ldots,c_{k,n}^{RBF},c_{k,1}^{p},\ldots,c_{k,M}^{p}\in\mathbb{R}$ are chosen to satisfy the interpolation conditions $s(\mathbf{x}_{k,j})=f(\mathbf{x}_{k,j})$, $j=1,2,\ldots,n$, along with constraints $\sum_{j=1}^{n}c_{k,j}^{RBF}\pi_{l}(\mathbf{x}_{k,j})=0$, for $l=1,2,\ldots,M$.

2.3 Step 3: Integrate the interpolant of the integrand

By integrating the interpolant, the approximation of the integral of $f$ is reduced to

[TABLE]

for $k\in\mathcal{K}_{I}$ and

[TABLE]

for $k\in\mathcal{K}_{S}$. A simple derivation can be carried out to show that the weights can be found by solving the linear system $A_{k}\mathbf{w}_{k}=\mathbf{I}_{k}$ with $(n+M)\times(n+M)$ matrix

[TABLE]

The $n\times n$ submatrix $\Phi_{k}$ is made up of the RBFs evaluated at each point in $\mathcal{N}_{k}$, that is

[TABLE]

Likewise the $n\times M$ matrix $P_{k,il}$ consists of the polynomial basis evaluated at each point in $\mathcal{N}_{k}$ so that

[TABLE]

If $k\in\mathcal{K}_{I}$, the right hand side, $\mathbf{I}_{k}$, includes integrals of the basis functions over $t_{k}$ only. That is,

[TABLE]

The integrals of the trivariate polynomial terms can be evaluated exactly via, for instance, the Divergence Theorem or barycentric coordinates. For the RBFs, the integrals can be evaluated by further decomposing $t_{k}$ into four tetrahedra that share a common vertex. Summing the integrals of the RBFs over the four tetrahedra results in the integral over $t_{k}$. This process allows the volume integral over a tetrahedron to be reduced to four integrals in a single dimension. Section 2.3.1 explains this process.

On the other hand, for $k\in\mathcal{K}_{S}$

[TABLE]

The integrals over $t_{k}$ are evaluated using the methods described in the previous paragraph while the integrals over $s_{k}$ are approximated using a scheme discussed in section 2.3.2.

2.3.1 Integrals of RBFs Over Tetrahedra

Suppose that the tetrahedron $t_{k}$ has vertices $\mathbf{a}_{k}$, $\mathbf{b}_{k}$, $\mathbf{c}_{k}$ and $\mathbf{d}_{k}$, all points in $\mathbb{R}^{3}$. Let $\mathbf{x}_{k,j}$ be some point in $\mathbb{R}^{3}$, which will be common to the four tetrahedra that will be integrated over to obtain the integral over $t_{k}$. Although what follows applies for any point in $\mathbb{R}^{3}$, the point $\mathbf{x}_{k,j}$ is an interpolation node from the set $\mathcal{N}_{k}$ in this context. A unit length normal vector to the side of $t_{k}$ with vertices $\mathbf{a}_{k}$, $\mathbf{b}_{k}$ and $\mathbf{c}_{k}$ is defined by

[TABLE]

and when defining a normal vector in what follows, the order of the vertices matters and should be taken as the order shown in this definition. Further, let $\mathbf{e}_{k,j}$, $\mathbf{f}_{k,j}$, $\mathbf{g}_{k,j}$ and $\mathbf{h}_{k,j}$ be the orthogonal projections of $\mathbf{x}_{k,j}$ onto the sides of $t_{k}$ with vertices $\mathbf{a}_{k}$, $\mathbf{b}_{k}$ and $\mathbf{c}_{k}$; $\mathbf{a}_{k}$, $\mathbf{d}_{k}$ and $\mathbf{b}_{k}$; $\mathbf{a}_{k}$, $\mathbf{c}_{k}$ and $\mathbf{d}_{k}$; and $\mathbf{b}_{k}$, $\mathbf{d}_{k}$ and $\mathbf{c}_{k}$, respectively. For instance,

[TABLE]

Then by applying the divergence theorem it can be shown that

[TABLE]

This expression for the integral over $t_{k}$ contains integrals over the four tetrahedra that share $\mathbf{x}_{k,j}$ as a common vertex. Consider

[TABLE]

since the remaining integrals over the tetrahedra are analogous. The integrand is radially symmetric about $\mathbf{x}_{k,j}$ and depends only on the distance from $\mathbf{x}_{k,j}$ suggesting the change of variables

[TABLE]

In this change of variables, triangles similar to the side of $t_{k}$ with vertices $\mathbf{a}_{k}$, $\mathbf{b}_{k}$ and $\mathbf{c}_{k}$ are parameterized using barycentric coordinates and scaled by the nonnegative parameter $\sigma$, which also accounts for the distance of the triangle from the point $\mathbf{x}_{k,j}$. Under this change of variables the integral becomes

[TABLE]

where $V_{k,j}=\left|(\mathbf{a}_{k}-\mathbf{x}_{k,j})\cdot\left[(\mathbf{b}_{k}-\mathbf{x}_{k,j})\times(\mathbf{c}_{k}-\mathbf{x}_{k,j})\right]\right|$ is six times the volume of the tetrahedron $t_{\mathbf{x}_{k,j}\mathbf{a}_{k}\mathbf{b}_{k}\mathbf{c}_{k}}$. This volume appears in the Jacobian determinant from the change of variables. Also,

[TABLE]

is the value of $\sigma$ corresponding to the side of $t_{k}$ with vertices $\mathbf{a}_{k}$, $\mathbf{b}_{k}$ and $\mathbf{c}_{k}$.

Now, in the case of $\phi(r)=r^{2p+1}$, $p=0,1,2,\ldots$, the iterated integrals in $\sigma$ and then $\lambda_{2}$ can be computed in closed form. However, exploring the integration over $\lambda_{1}$ in Mathematica indicates the cost of a closed form expression for the integral is computationally too expensive, so the proposed algorithm uses standard pseudospectral methods for evaluating the integrals over $\lambda_{1}$.

2.3.2 Integrals Over Slivers of Volume at the Surface

When assigning a sliver of volume to a particular tetrahedron, $t_{k}$, care must be taken so that there are no gaps or overlaps between adjacent slivers. Let $\tau_{k,i}$, $i=1,2,3,4$, be the triangular faces of $t_{k}$. At least one of these faces has all three vertices on the surface of the sphere. In most cases, this will be only one face of $t_{k}$ (particularly when the volume is well resolved by small enough tetrahedra), call it $\tau_{k,*}$. It turns out that if the three edges of $\tau_{k,*}$ are projected radially from the center of the sphere to the surface, gaps and overlaps will be prevented. For each edge of $\tau_{k,*}$ the area between the arc on the sphere surface and the edge of the triangle forms a side of the sliver of volume. The boundary of the sliver volume is formed by all three of these sides, the spherical triangle on the surface of the sphere between the three sides, and the triangle $\tau_{k,*}$. Figure 1 illustrates one of these volumes.

Assigning the slivers of volume in this way provides for a transformation of the coordinates of the sliver which allows the integral over the volume to be written as an iterated integral over a triangular area and a parameter, $\sigma$, which relates to the projection from the origin. Consider any point $\mathbf{x}$ in the volume. The vector $\mathbf{x}$ intersects the plane containing $\tau_{k,*}$ at a point $\mathbf{x}^{\prime}$ that is inside the triangle $\tau_{k,*}$. All of the points inside $\tau_{k,*}$ can be parameterized by, for instance,

[TABLE]

where $\mathbf{a}_{k}$, $\mathbf{b}_{k}$ and $\mathbf{c}_{k}$ are now representing the vertices of $t_{k,*}$. With this parameterization of $t_{k,*}$ any point $\mathbf{x}$ in the volume of the sliver can be represented as

[TABLE]

where $0\leq\sigma\leq\rho-\left\lVert\mathbf{x}^{\prime}(\lambda,\mu)\right\rVert_{2}$. Here $\sigma$ measures the distance from $\mathbf{x}^{\prime}$ to $\mathbf{x}$. With this parameterization, for instance,

[TABLE]

where $J$ is the Jacobian determinant of $\mathbf{x}$ with respect to $\sigma$, $\lambda$ and $\mu$. The integrals in $\sigma$, $\lambda$, and $\mu$ can be easily treated with any number of quadrature methods over an interval.

2.4 Step 4: Combine weights from the subdomains

Summing over all $k\in\{1,2,\ldots,K\}$ leads to the approximation of the volume integral over $\Omega$

[TABLE]

Let $\mathcal{K}_{i}$, $i=1,2,\ldots,N$, be the set of all pairs $(k,j)$ such that $\mathbf{x}_{k,j}\mapsto\mathbf{x}_{i}$. Then the volume integral over $\Omega$ can be rewritten as

[TABLE]

3 Test Examples

To demonstrate the performance of the method described herein, the algorithm will be applied to four different test integrands featuring varying degrees of smoothness. The fourth of these test integrands includes extremely localized feature. Weights are computed on quasi-uniformly spaced nodes, pseudo-randomly generated nodes, and a clustered node set with increased density near the localized feature of the fourth test integrand. This clustered node set is used demonstrate the performance of the algorithm under localized node refinement. In all of these tests, the radius of the ball is fixed to $\rho=\frac{6^{\frac{1}{3}}}{2\pi^{\frac{1}{3}}}$ so that the volume of the ball is equal to one.

3.1 Node Sets

In the cases of quasi-uniformly spaced nodes and the clustered node set, quadrature nodes were generated using a modification of the algorithm presented in [14]. The pseudo-randomly spaced node sets were generated by first drawing a set of points from the two-dimensional Halton sequence and mapping the set to the surface of the ball. Then points were drawn from the three-dimensional Halton sequence, mapping the set to the domain $(x,y,z)\in[-\rho,\rho]\times[-\rho,\rho]\times[-\rho,\rho]$ and keeping only points satisfying $x^{2}+y^{2}+z^{2}\leq\rho-\frac{h}{10}$, where $h$ is prescribed to be the average spacing between the nodes on the surface. Examples of the node sets are displayed in figure 2.

3.2 Performance on Test Integrands

The algorithm was applied to four test integrands. When generating quadrature weights, in all cases the radial basis function was $\phi(r)=r^{3}$ and the number of nearest neighbors, $n=(m+1)(m+2)(m+3)$, was based on the trivariate polynomial order $m$. Some computational experiments in, for instance, [17, 3, 2] indicated that in the presence of boundaries the number of nearest neighbors must be large enough to overcome effects like Runge phenomenon. The examples given in [17] indicated that the boundary errors were most prominent when nodes were (exactly) uniformly spaced. Therefore, to determine how many nearest neighbors should be included to overcome boundary errors, the algorithm here was modified to compute volume integrals over a cube. When considering a cube the entire volume can be decomposed by tetrahedra, so the algorithm need not consider slivers of volume near the surface. Figure 3 illustrates the absolute error when integrating $f(x,y,z)=\frac{1}{1+((x-x_{s})^{2}+(y-y_{s})^{2}+(z-z_{s})^{2})}$, with $x_{s}=0.234841098236337$, $y_{s}=0.048716273957102$ and $z_{s}=0.214415743035283$, over the volume of the unit cube centered at the origin for various choices of $n$ and $m$. The matrix $A_{k}$ is singular for choices of $n$ below $n=\frac{(m+1)(m+2)(m+3)}{6}$, this is indicated by the lower dashed curve in the figure. Further, for each case of $m=0,1,2,\ldots,7$ it is clear that choices of $n$ below $n=(m+1)(m+2)(m+3)$ can lead to large errors.

Performing a similar experiment for the unit ball, figure 4 illustrates that even the relaxation from nodes that are exactly uniformly spaced to those that are quasi-uniformly spaces can allow for $n$ as low as $(m+1)(m+2)(m+3)$. However, all results shown here utilize $n=(m+1)(m+2)(m+3)$.

The first of the test integrands is a degree 30 trivariate polynomial. That is, let

[TABLE]

The coefficients of the polynomial are available in a Matlab file at [15]. The exact value of the integral over the ball is

[TABLE]

with $\Gamma$ the gamma function, and for the set of coefficients used here the expression evaluates to 3.792079311949332. Figure 5 displays convergence of the approximate integral to the exact value at an order better than O$\left(N^{-\frac{m}{3}}\right)$, where $m$ corresponds to the order of the trivariate polynomial terms used in the approximation. If $h$ refers to a typical node separation distance, this corresponds to a convergence order of better than O$\left(h^{m}\right)$, especially in the case of quasi-uniformly spaced nodes. The theory in [1] explains that if the multivariate polynomial basis up to degree $m$ is included in the process of RBF interpolation, then all of the terms in the Taylor series up to degree $m$ will be handled exactly for the function being interpolated. The remaining terms in the Taylor series are then approximated by the RBF basis that was included. This is leading to a convergence order of at least O$(h^{m})$.

The second test integrand is the Gaussian

[TABLE]

where

[TABLE]

is a randomly chosen shift of the center of the Gaussian from the origin. In order to have an accurate value to compare to, the volume integral was first approximated by evaluating

[TABLE]

using Matlab’s integral3 command with the absolute and relative tolerances both set to ten times machine precision. The resulting approximation of the integral for $f_{2}$ is 0.161965667295343. Figure 6 illustrates the error in the integral of $f_{2}$ over the ball when compared to the result from Matlab after rotating the integrand randomly 100 times. It is clear again that the order of the error is most dependent on the degree of the polynomials used in the interpolation.

The third test integrand is

[TABLE]

with sign the signum function. This function is discontinuous at the plane $z=0$, so any method based on a continuous approximation of the integrand across the discontinuity should not be expected to achieve better than O$\left(N^{-\frac{1}{3}}\right)$ (i.e. O$(h)$) error. Figure 7 illustrates this for the present method.

3.3 Performance When Utilizing Clustered Node Sets

To illustrate further utility of the proposed method, the algorithm was also applied to a test integrand featuring a steep and localized gradient. To capture the rapid change in the integrand node sets were generated that feature more densely clustered nodes near the local feature. The test integrand was

[TABLE]

which has a steep gradient near the origin. Figure 8 illustrates that in cases where quasi-uniformly spaced or pseudo-randomly spaced node sets cannot capture the changes in the integrand, weights generated for node sets with clustering near features improve the approximation under refinement.

3.4 Computational Expense

Just like the algorithms presented in [16, 18, 17], the ability to consider each tetrahedron individually allows the time to compute a set of quadrature weights and the use of memory both to scale like $O(N)$. Figure 9 illustrates the time to compute the set of quadrature weights on $N$ nodes for various choice of the polynomial order, $m$. Since the choice of $m$ affects the sizes of the systems of linear equations that need to be solved at each iteration, the figure shows an increase in the computational cost as $m$ increases.

Further, except for the identification of nearest neighbors in order to construct the local weight set for each tetrahedron/sliver of volume and for the combination of weights in step 4 the algorithm is pleasingly parallel. The parallelization tests in [16] illustrate that these two steps do not have a significant impact on the scalability of the algorithm with the number of cores when considering evaluating surface integrals, and the same is true for the volume integration algorithm described herein.

4 Conclusions

This study has supplemented the previous RBF-FD based approach for evaluating definite integrals [16, 18, 17] with an extension to integrals over volumes. The computational tests illustrate an algorithm that can achieve at least $O(h^{m})$ accuracy, with $h$ the typical node separation distance and $m$ the order of trivariate polynomial basis functions included in the approximation. On a set of $N$ nodes in the ball, the computational cost is only $O(N)$ and the algorithm is pleasingly parallel. A key feature of the algorithm is that it is able to compute quadrature weights on even irregularly spaced or clustered node sets.