Using Cubature in python

Question

I am trying to calculate a seemingly simple integral. Here is the MMA code:

In[1]:= rect[x_] = If[x <= 0.5, 1, 0];

In[2]:= N[Integrate[Integrate[rect[(x^2 + y^2)^.5], {x, 0, 100.}], {y, 0., 100.}]]

Out[3]= 0.19635

However, when I do that with the cubature package in Python, I cannot get anything reasonable for the integration limits over 13, and even for (-13,13) it lies. It really has to be something like +/- 10.

Here is the python code:

import time
import numpy as np
from cubature import cubature

start = time.time()


def rect(r):
    return np.where(abs(r) <= 0.5, 1, 0)


def Wp(r1, r2):
    return rect((r1 ** 2 + r2 ** 2) ** 0.5)


def integrand_rectangle(x_array):
    return Wp(x_array[:, 0], x_array[:, 1])


# boundaries
a = 0.0
b_x = 100.0
xmin_t = [a, a]
xmax_t = [b_x, b_x]

val_t, err_t = cubature(integrand_rectangle, 2, 1, xmin_t, xmax_t, vectorized=True)
print(val_t)

end = time.time()

print(f"Runtime of the program is {end - start}")

I have tried to play with the number of interactions, relative and absolute error tolerances, nothing seems to cure this. Is there a solution? Why can't this well-tested package perform such a simple integration? Am I doing it wrong?

I have tried original C-cubature and it works perfect and fast! It is a problem of the python -- C interface! It is a bug!

Take a look at [quadpy](https://github.com/nschloe/quadpy/). — Nico Schlömer, Sep 22 '20 at 13:59
It looks like you want to integrate the area of a disk of radius `0.5`. You wouldn't do that using quadrature. — Nico Schlömer, Sep 22 '20 at 14:04
Original cubature package in C does it accurate and fast. It is only cubature wrapped for python lies. I am not good enough in python to spot where exactly the interface is corrupted. Also, I tried quadpy, nquad and mcint. First two are inaccurate and slow. Third one produces garbage. What is a good package to use then? — MsTais, Sep 22 '20 at 14:34
What problem are you trying to solve? The area of a something. Quadrature tools are not suitable for the task. — Nico Schlömer, Sep 22 '20 at 14:51
Can you try the C-cubature? It does it correctly. Suitable or not, the original C-package does it. The one with python - C interface fails. So this is not cubature itself. This is the interface. — MsTais, Sep 22 '20 at 14:59
To rephrase: the original C-package and the C-package-based python package produce different results, where C-package is right and python package is wrong. Python package uses C-package + interface. Hence, the problem is NOT the package, but the interface. I hope it makes better sense now. — MsTais, Sep 22 '20 at 15:02
On top of that "not-suitalbe"... Why, what are the limitations? Cubature is the adaptive integration procedure. Adaptive algorithms are one of the only options for this integral, since it is so delta-function-like. It is worth it to explain what you mean by that as well... — MsTais, Sep 22 '20 at 15:05

score 1 · Accepted Answer · answered Sep 22 '20 at 15:17

It looks like you're trying to integrate a constant function over a domain (a disk of radius 1/2 in your case) to find out the area/volume of the domain. Although the problem can be formulated as an integral,

one would not typically do that with quadrature. The reason for this is that most quadrature methods (including the one you use) are Gaussian quadrature methods, meaning that they will produce approximate results. The approximation is very good for smooth functions (i.e., polynomials of low degree). The function you have is not smooth at all, and even discontinuous at the boundary of the domain. That's why quadrature results will be inaccurate.

A more appropriate approach would be to try and triangulate the domain, using one of the many available tools.

After meshing, you can add the volumes of the triangles to get the domain volume.

Using Cubature in python

1 Answers1