您如何像Mathematica一样执行这种不合适的积分？

您如何像Mathematica一样执行这种不合适的积分？

一个非常有趣的问题。

首先要注意的是

from numpy import expdef f(x):    return exp(-x)def g(x):    c = 0.9    return c * x**(c - 1) * exp(-x ** c)def integrand(x):    return f(x) * log(f(x) / g(x))

在0处具有 可积分 的奇点，并且[0，infty]上的积分可以通过 分析 求值。经过一些 *** 作后，您会发现

import numpyimport scipy.specialc = 0.9# euler_mascheroni constantgamma = 0.57721566490153286060val = scipy.special.gamma(c + 1) - 1 - numpy.log(c) + (c - 1) * gammaprint(val)0.0094047810750603

wolfram-
alpha
正确地将其值赋予许多数字。为了用数值方法重现这一点，一个好的尝试总是总是tanh-
sinh正交
（例如，来自我的一个项目Quadpy）。以某个较大的值切断该域，无论如何该函数几乎都为0，然后：

from numpy import exp, logimport quadpydef f(x):    return exp(-x)def g(x):    c = 0.9    return c * x**(c - 1) * exp(-x ** c)def integrand(x):    return f(x) * log(f(x) / g(x))val, err = quadpy.tanh_sinh(integrand, 0.0, 100.0, 1.0e-8)print(val)0.009404781075063085

---

现在，对于其他一些事情，也许令人惊讶的是，它做得 不太 好。

当看到这种类型的整数时

exp(-x) * f(x)

，首先想到的是高斯-
拉格瑞正交
。例如，使用Quadpy（我的一个项目）：

import numpyimport quadpyc = 0.9def f(x):    return numpy.exp(-x)def g(x):    return c * x ** (c - 1) * numpy.exp(-x ** c)scheme = quadpy.e1r.gauss_laguerre(100)val = scheme.integrate(lambda x: numpy.log(f(x) / g(x)))print(val[0])

这给

0.010039543105755215

尽管我们使用了100个积分点，但对于实际值而言却是一个令人惊讶的差值。这是由于以下事实：多项式，尤其是项

log(x)

和不能很好地逼近被积数

x **c

：

import numpyfrom numpy import exp, log, onesfrom scipy.special import gammaimport quadpyc = 0.9def integrand(x):    return exp(-x) * (-x - log(c) - (c - 1) * log(x) - (-x ** c))scheme = quadpy.e1r.gauss_laguerre(200)val = scheme.integrate(lambda x: -x - log(c) - (c - 1) * log(x) - (-x ** c))[0]vals = numpy.array([    - scheme.integrate(lambda x: x)[0],    -log(c) * scheme.integrate(lambda x: ones(x.shape))[0],    -(c - 1) * scheme.integrate(lambda x: log(x))[0],    scheme.integrate(lambda x: x ** c)[0]])euler_mascheroni = 0.57721566490153286060exact = numpy.array([    -1.0,    -log(c),    euler_mascheroni * (c-1),    gamma(c + 1)])print("approximation, exact, diff:")print(numpy.column_stack([vals, exact, abs(vals - exact)]))print()print("sum:")print(sum(vals))approximation, exact, diff:[[-1.00000000e+00 -1.00000000e+00  8.88178420e-16] [ 1.05360516e-01  1.05360516e-01  6.93889390e-17] [-5.70908293e-02 -5.77215665e-02  6.30737142e-04] [ 9.61769857e-01  9.61765832e-01  4.02488825e-06]]sum:0.010039543105755278