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 个数量级。提前谢谢您。

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

Answer 1

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

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

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

我编写了一个名为

Answer 2

的 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