10

as the title suggests I'm trying to compute the integral of a compactly supported function (Wendland's quintic polynomial) on a triangle. Notice, that the center of the function is somewhere in 3-D space. I integrate this function on an arbitrary, but small triangle ($area < \frac{(radius/4)^2}{2}$). I am currently using the integration described by Dunavant, 1985 (p=19).

It seems however, that these quadrature rules are not suited towards compactly supported problems. This is supported by the fact that when I integrate $f(r) = [r\leq1]$ (so a function that is 1 inside the circle of radius 1) on a plane which is discretized using triangles, my (normalized) results are between 1.001 and 0.897.

So my question is, does a specialized quadrature rule exist for this kind of problem? Would a lower order composite integration rule work better?

Unfortunately this routine is really critical in my code so precision is crucial. On the other hand I need to do this integration "a couple of times" for a single time-step so computational expense should not be too high. Parallelization is not an issue as I will execute the integration itself in serial.

Thanks in advance for your answers.

EDIT: Wendland's quintic polynomial is given by $W(q) = [q\leq2]\frac{\alpha}{h^3}(1-\frac{q}{2})^4(2q+1)$ with $\alpha = \frac{21}{16\pi}$ and $q=\frac{\|r-r_0\|}{h}$ with $r_0$ being an arbitrary vector in $\mathbb{R}^3$

EDIT2: If $\Delta$ is the two-dimensional triangle then I want to calculate $\int_\Delta \omega(r) dr$ with $\omega(r) = W(\frac{\|r-r_0\|}{h})$. So $q$ in $W$ will never be smaller than 0. Note that the integral is a surface integral over a 2-D surface in $\mathbb{R}^3$

EDIT3: I have an analytical solution for the 1-D (line) problem. Calculating one for 2-D (triangle) might be possible as well.

Azrael3000
  • 307
  • 1
  • 8

3 Answers3

4

Since the function is smooth within $q \leq 2$, but not of fixed degree (in the plane, that is), I would suggest using a simple adaptive scheme, e.g. the Trapezoidal Rule with Romberg's method, in both dimensions.

That is, if your triangle is defined by the vertices $x$, $y$ and $z \in \mathbb R^3$, and you have a routine romb(f,a,b) which integrates f along the line from a to b, you could do the following (in Matlab notation):

int = romb( @(xi) romb( W , xi , y+(z-y)*(xi-x)./(z-x) ) , x , z );

In romb, don't use a fixed number of points, but keep growing the table until the difference between two successive diagonals is below your required tolerance. Since your function is smooth, this should be a good error estimate.

If parts of the triangle are outside of the domain of $W(q)$, you could try adjusting the limits of integration in the above code accordingly.

This may not be the most computationally efficient way of solving your problem, but the adaptivity will give you much more robustness than a fixed-degree rule will.

Pedro
  • 9,573
  • 1
  • 36
  • 45
  • The function is smmoth everywhere except for $q=0$. The neighborhood of this point is causing the trouble. – Arnold Neumaier Mar 22 '12 at 14:15
  • Ah decomposing into two 1-D problems, not a bad idea at all. Because there is one thing that I haven't told you. I have an analytical solution in 1-D so I can replace the interior romb by an analytical function. I'll give that a shot +1 already – Azrael3000 Mar 22 '12 at 14:16
  • @ArnoldNeumaier, I'm sorry, I don't see how that is possible. Could you explain? – Pedro Mar 22 '12 at 14:45
  • smooth as a function of $q$, but $q$ is a nonsmooth function of $r$, and integration is over $r$, as far as I understood the question. The composite function is thus a nonsmooth function of $r$. – Arnold Neumaier Mar 22 '12 at 16:50
  • @ArnoldNeumaier, ok, I see, but this will only be a problem if $r$ is in the two-dimensional triangle. As soon as it is not, $q>0$. – Pedro Mar 22 '12 at 16:56
  • $q>0$ but small still behaves numerically poorly. Any numerical method for any problem is unstable close to a singularity, not only at the singularity. See my updated answer for details. – Arnold Neumaier Mar 22 '12 at 17:03
  • 1
    @Pedro I implemented it and it works like a charm. We actually also found an analytical solution today. But this is only for a special case which can be used to reconstruct the general one. That means we need to do some domain decomposition. Since the Romberg converges in about 4 steps I think that because of this it will be faster than using the analytical formula. And according to Wikipedia we can do still better than Romberg when using rational polynomials. You'll find your name in the acknowledgements of my next paper :) Cheers. – Azrael3000 Mar 30 '12 at 22:02
  • @Azrael3000: Glad to hear it worked! Rational Polynomials for convergence accelerations sounds like the $\varepsilon$-Algorithm, which will need at least three steps to do an extrapolation (if I remember correctly), and it can be quite unstable. Thanks for the Ack though! Would be cool if you could mail me a link to the paper once it's out there :) – Pedro Mar 30 '12 at 23:45
  • We will try to submit something in summer. We'll see how long it takes until it goes through. But I'll send you a copy once it's out. – Azrael3000 Mar 31 '12 at 09:51
2

Integration rules assume that the function is locally well approximated by a low degree polynomial. Your problem has nothing to do with compact support. Compactly supported radial basis functions are smooth at the support boundary, and quadrature rules up to the order of smoothness can be used without problems. (Higher order rules do not help; thus you should probably not use a rule that integrates degree 5 polynomials exactly.)

In your case, the inaccuracy comes from the fact that the assumption of good polynomial approximability fails in your case for triangles near $r_0$, even when they do not contain $r_0$.

$W$ is smooth as a function of $q$, but $q$ is a nonsmooth function of $r$, with a gradient that becomes infinite in the limit $r\rightarrow r_0$. Integration is over $r$, and the composite function is a nonsmooth function of $r$.

If the triangle does not contain $r_0$, the function is $C^inf$ but this doesn't help as the higher derivative grow very quickly close to $r_0$, and the error of a high order method is proportional to a high order derivative, hence very large!

The simple remedy is to split each triangle T into a number N_T of subtriangles. You can take $N_T=1$ far away from $r_0$, and $N_T\gg 1$ close to $r_0$. You can figure out offline how big $N_T$ must be for triangles of a given diameter and distance from $r_0$ to reach a desired accuracy. Moreover, you should only use low order formulas close to $r_0$.

As you integrate over a triangle, but $r_0$ is 3-dimensional, the triangle is apparently in $R^3$.

A faster remedy would therefore tabulate the integral for $r_0=0$ as a function of the triangle coordinates (normalized by rotating it into a 2-dimensional $xy$-plane such that one vertex lies on the $x$-axis, and reflecting it such that a second vertex lies above it). This tabulation must be sufficiently detailed to make a linear or quadratic interpolation accurate enough. But you can use the slow method outlined first to create this table.

Another way to get rid of the problem is to use a compactly supported radial basis function that is a polynomial in $q^2$ rather than $q$. This is smooth everywhere, and easy to integrate.

Arnold Neumaier
  • 11,318
  • 20
  • 47
  • I think there is a small missunderstanding. I updated the description of my question. As a matter of fact, in the integral $q$ can never be smaller than 0. And $r_0$ is not necessarily contained in the triangle. – Azrael3000 Mar 22 '12 at 14:35
  • Your new addition doesn't make sense to me. If $r_0\in R^3$ then so must be $r$. Or do you integrate over a 2D triangle in $R^3$? - I didn't assume that $r_0$ is in the triangle. I just added in a moment more detail to my answer. – Arnold Neumaier Mar 22 '12 at 16:52
  • Yes, it is correct that I integrate over a 2D triangle in $R^3$. – Azrael3000 Mar 22 '12 at 19:06
2

For a good overview of cubature rules, see "R. Cools, An Encyclopaedia of Cubature Formulas J. Complexity, 19: 445-453, 2003". Using a fixed rule, can give you the advantage that some rules integrate polynomials exactly (as Gaussian quadrature does in one-dimension).

Cools is also one of the main authors of CUBPACK, a software package for numerical cubature.

GertVdE
  • 6,179
  • 1
  • 21
  • 36