在 python 中使用 Cubature

Question

我正在尝试计算一个看似简单的积分。这是 MMA 代码：

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

但是，当我使用 Python 中的 cubature 包执行此操作时，对于超过 13 的积分限制，甚至 (-13,13)，我都无法得到任何合理的信息。它确实必须是 +/- 10 之类的东西。

这里是 python 代码：

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}")

我尝试过交互次数、相对和绝对误差容限，似乎没有什么可以解决这个问题。有解决办法吗？为什么这个经过良好测试的包不能执行如此简单的集成？我做错了吗？

我已经尝试过原始的 C-cubature，它完美且快速！ python——C接口的问题！这是一个错误！

Answer 1

看起来您正试图在一个域（在您的例子中是一个半径为 1/2 的圆盘）上积分一个常量函数来找出域的 area/volume。虽然问题可以表述为一个积分，

人们通常不会用正交来做到这一点。这样做的原因是大多数正交方法（包括您使用的方法）都是 Gaussian 正交方法，这意味着它们将产生近似结果。该近似值非常适用于平滑函数（即低阶多项式）。你的函数一点也不平滑，甚至在域的边界处不连续。这就是正交结果不准确的原因。

更合适的方法是使用 one of the many available tools.

尝试对域进行三角测量

网格化后，可以将三角形的体积相加得到域体积。

在 python 中使用 Cubature

Using Cubature in python

python

numerical-methods