python - 计算积分的有效方法？

Question

我知道如何使用蒙特卡洛计算积分，但我想知道是否可以使用梯形法则结合 numpy 提高效率以获得相同的积分，我不确定哪个最快或后者可能吗？

例如集成 e**-x**2 > y 我可以像这样使用蒙特卡洛方法:

import numpy as np
import matplotlib.pyplot as plt

X = np.random.rand(500000,2)
X[:,0] = X[:,0]*4-2
J = np.where(X[:,1] < np.exp(-X[:,0]**2))[0]
Xp = X[:2000]
Ip = [i for i in range(len(Xp)) if i in J]
Inp = [i for i in range(len(Xp)) if i not in J]
plt.plot(Xp[Ip,0],Xp[Ip,1], 'bd', Xp[Inp,0],Xp[Inp,1], 'rd')
plt.show()

这很容易计算:

print len(J) / 500000.0 * 4

给出:

1.767784

在这种情况下很容易，但如果没有指定间隔，如 [a,b] , n，我想创建一个函数，那么我认为上面的方法不是真的高效，至少我是这么认为的。

所以，我的问题是，我能否将一个连续函数(例如 cos(x)/x)集成到一个函数中，例如 [a,b]梯形规则？

它比我在这里使用的方法更好吗？

欢迎大家提出建议。

Answer 1

只需使用scipy.integrate.quad:

from scipy import integrate
from np import inf
from math import exp, sqrt, pi

res, errEstimate = integrate.quad(lambda x: exp(-x**2), -inf, +inf)

print(res)       #output: 1.7724538509055159
print(sqrt(pi))  #output: 1.7724538509055159

最后一行简单地检查评估的积分确实是 Pi 的平方根(它是 Gaussian integral )。