scipy 四边形似乎提供不准确的积分值

Question

令我感到困惑的是，scipy 中的

integrate.quad

即使对于简单的积分也会返回不准确的结果。下面，我将给出一个简单的例子，可以很容易地与分析计算进行交叉检查。

但是，我想通过了解 scipy/numpy 中最好和最可靠的积分方法来了解该问题的一般解决方案。这是因为，我需要做的实际积分非常复杂，并且解析解和/或重新参数化是不可能的。

这是一个简单的例子：

from scipy.integrate import quad
xmax=1.0e23
xmin=0.5
def myFUN(x):
    return 1/(1+x)**(3.0)
print( quad( myFUN, xmin,xmax )[0] )

返回：

2.3171047961935996e-40

现在我们来分析一下：

def ana(x):
    return -0.5/(1+x)**2.0 
print( ana(xmax)-ana(xmin) )

一个人得到：

0.2222222222222222

请注意，数值输出与分析计算相差近 40 个数量级。

提前谢谢您。

Answer 1

我编写了一个名为 lintegrate 的 Python 集成包，它适用于函数的对数，您也可以使用它。安装 lintegrate 后：

pip install lintegrate

你可以这样做：

import numpy as np

from lintegrate import lcquad

def myfun(x, *args):
    """
    Define the natural logarithm of the function you want to integrate.
    """

    # natural log of 1 / (1 + x)**3
    return -3.0 * np.log(1 + x)


xmin = 0.5
xmax = 1e23

# return the natural log of the original function
lint = lcquad(myfun, xmin, xmax)[0]

# exponentiate to get your required answer
print(f"Integral = {np.exp(lint)}")
Integral = 0.2222222222800475

您可以在

lcquad

此处查看文档。它在内部使用 GSL cquad 函数。

Answer 2

想象一下你正在用梯形规则来解决这个问题。您将范围分成 1000 个间隔。如果 xmin = 0.5 且 xmax = 1e23，则第一个内部值为 1e20，函数的值约为 1e-60，或者，对于所有意图和目的，0。显然 scipy 的四边形例程将使用比梯形规则更高级的东西，但是后果是一样的。

因此，如果您选择使用以“近无穷大”作为参数之一的库例程，那么您将必须采取特殊措施（例如，将积分范围截断到可以将被忽略的部分限制在可容忍范围以下的范围）水平）。

如果您使用合理的积分范围，则该函数表现良好。

from scipy.integrate import quad
xmax=10.0
xmin=0.5

def myFUN(x):
    return 1/(1+x)**(3.0)

def indefinite_integral( x ):
    return -0.5 / ( 1 + x ) ** 2

print( "Using library routine: ", quad( myFUN, xmin,xmax )[0] )
print( "Using your maths knowledge: ", indefinite_integral( xmax ) - indefinite_integral( xmin ) )

输出：

Using library routine:  0.21808999081726355
Using your maths knowledge:  0.21808999081726355

编辑：

通常可以对积分进行变换，使积分范围为 O(1)。例如，在您的情况下，您可以替换

u = 1/(1+x)³

然后 f(x).dx 变换为 (-1/3)u^-1/3.du.

在下面的代码中，我回到了原来的积分限制，但也进行了此替换。 （我还翻转了积分限制以使被积函数为正，但这并不重要。）然后它会在整个范围内为您提供您想要的答案。

from scipy.integrate import quad
xmax=1e23
xmin=0.5

def myFUN( x ):
    return 1.0 / ( 1.0 + x ) ** 3.0

def indefinite_integral( x ):
    return -0.5 / ( 1 + x ) ** 2

print( "Using library routine: ", quad( myFUN, xmin,xmax )[0] )
print( "Using your maths knowledge: ", indefinite_integral( xmax ) - indefinite_integral( xmin ) )

######################

def altFUN( u ):
    return 1.0 / ( 3.0 * u ** ( 1.0 / 3.0 ) )

umin = myFUN( xmax )
umax = myFUN( xmin )
print( "Using transformed integral: ", quad( altFUN, umin, umax )[0] )

输出：

Using library routine:  2.3171047961935996e-40
Using your maths knowledge:  0.2222222222222222
Using transformed integral:  0.2222222222222222