高效计算积分的方法？

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]？

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

欢迎任何建议。

- JKM

---

嗯，不，将一维的、边界良好的、具有快速衰减导数的可微函数与蒙特卡罗方法集成确实是一个非常糟糕的想法。但你从哪里得到梯形法则也不是完全垃圾的想法呢？还有更好的数值积分方法。scipy.integrate有什么问题吗？ - Andrey Tyukin

@AndreyTyukin 我在编程方面还是个新手，我发现其他方法使用了太多新命令难以理解，所以我从最简单的开始。但如果你有任何特定的方法推荐，我会认真听取。 - JKM

"非常初学者"是相对的......一个普通的“初学者”能否实现蒙特卡罗方法计算积分并不明显。这也不是特定于编程的。 - Andrey Tyukin

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的平方根（它是高斯积分）。

- Andrey Tyukin

Answer 2

你也可以使用黎曼逼近法。下面的代码是Java语言实现的

package math;

import java.util.Optional;
import java.util.function.*;
import java.util.stream.IntStream;
import static java.lang.Math.*;

public class IntegralJava8
{
    public interface Riemann extends BiFunction<Function<Double, Double>, Integer, 
    BinaryOperator<Double>> { }

    public static void main(String args[])
    {
        int N=100000;
        Riemann s = (f, n) -> (a, b) -> 
        IntStream.range(0, n)
        .mapToDouble(i -> f.apply(a + i * ((b - a) / n)) * ((b - a) / n)).sum();
        Optional<Double> gaussIntegral = 
            Optional.of(s.apply(x -> exp(-pow(x, 2)), N).apply(-1000.0, 1000.0));
        gaussIntegral.ifPresent(System.out::println);
    }
}

在上述类中，它将计算从负无穷到正无穷的高斯积分，其值等于圆周率的平方根（1.772）。

- Jose Luis Soto Posada