1

我有一个问题scipy.quad。简而言之，我有一组非常长且复杂的嵌套函数和积分，其中包括必须在特定范围 10^2 < x < 10^20 内积分的递减函数的积分。

为了简单地演示这个问题，考虑使用 y=x^(-2) 在这些值之间的积分numpy.quad：

    import numpy as np
    from scipy.integrate import quad

    def func(x,z):
        "decreasing function to test"
        #print x,x**-2.
        return x**(-2.)

    #Test quad:
    print 'Quad', quad(lambda x: func(x,0.),1e2,1e20)[0]

这个积分的真正答案是 0.01，但是 quad 返回一个非常小的值，并且通过取消注释函数中的打印行，您可以看到它只对最大的 x 值进行积分（或者由于x 轴）。我需要找到一些方法来解决这个问题。

我知道您可以通过使用其他方法（例如辛普森一家或梯形规则）得到正确答案：

    from scipy.integrate import simps
    from numpy import trapz

    xarray=np.logspace(2,20,1000)

    print 'Simpson logspace x',simps(func(xarray,0.),xarray)
    print 'Trapezoid logspace x',trapz(func(xarray,0.),xarray)

但是，这些包括将数组传递给函数，这在我的实际代码中是不可能的。因此，我必须包含一个 for 循环来生成要与之集成的 y 数组，这会不可接受地减慢整个程序的速度。

有什么技巧可以让 quad 用于这种 x 范围，或者像 quad 一样工作的东西，我可以用它来代替吗？

4

0 回答 0