在python中使用Cubature

问题描述

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

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

In[2]:= N[Integrate[Integrate[rect[(x^2 + y^2)^.5],{x,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,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,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接口的问题！这是一个错误！

解决方法

您似乎正在尝试在某个域（您的情况下为半径为1/2的磁盘）上集成常数函数，以找出该域的面积/体积。尽管可以将问题表达为一个整体，

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

更合适的方法是使用one of the many available tools尝试对域进行三角测量。

网格划分后，您可以添加三角形的体积以获得域体积。