与 Python 的主要价值集成

Question

我正在尝试将一个函数与python集成。它有一个奇点，并且使用 scipy.integrate 我无法获得令人满意的结果（我还得到了 IntegrationWarning）。

这是我迄今为止的代码。这是我第一次处理这样的功能，我有点失落。

from scipy.integrate import quad
import numpy as np

def integrand(x):
    return -(x * np.sin(2*x)) / (x * np.cos(x) + np.sin(x))

def cauchy_principal_value(f, a, b, singular_point):
    integral = quad(f, a, b, points=[a, singular_point, b],limit=1000)
    return integral

result = cauchy_principal_value(integrand, 0, np.pi, 2.0287665755744362)

print("Principal Value: ", result[0])

有人可以帮助我吗？

Answer 1

一般来说，积分器例程在处理奇点时需要一些帮助。

在这里，如果您知道您想要一个主值，请使用

quad

告诉它

weight="cauchy"

（查找

quad

的文档）。这是一个快速而肮脏的演示：

In [32]: x0 = brentq(lambda x: x * np.cos(x) + np.sin(x), 1e-2, np.pi)

In [33]: quad(lambda x: -(x * np.sin(2*x)) / (x * np.cos(x) + np.sin(x) * (x - x0)), 1e-4, np.pi, weight='cauchy', wvar=x0)
Out[33]: (-4.14269240855162, 1.4398780224416236e-08)

• 积分器需要奇点的位置，所以用寻根器找到它（下面的

  brentq

  ）
• 与

  weight="cauchy"

  被积函数解释为

  f(x)

  中的

  f(x) / (x - c)

  。因此，围绕

  x0

  计算展开式并重新定义被积函数。下面我乘以

  x - x0

  来代替。
• x->0 处存在 0/0，因此理想情况下将

  integrand

  函数重新定义为零。下面我只是退出

  x=0

  。

Answer 2

scipy.integrate.quad 函数是通用积分器，在处理奇点或高度振荡函数时可能会遇到困难。在您的情况下，您正在尝试将函数与奇点集成，这可能会导致 IntegrationWarning 或不准确的结果。

要处理被积函数中的奇点，您可以使用柯西主值方法。这是代码的更新版本：

from scipy.integrate import quad
import numpy as np

def integrand(x):
    return -(x * np.sin(2*x)) / (x * np.cos(x) + np.sin(x))

def cauchy_principal_value(f, a, b, singular_point, **kwargs):
    def wrapped(x):
        if x == singular_point:
            return np.nan
        return f(x)
    
    integral = quad(wrapped, a, b, **kwargs)
    return integral

result = cauchy_principal_value(integrand, 0, np.pi, 2.0287665755744362)

print("Principal Value:", result[0])

在更新的代码中，引入了一个包装函数wrapped来处理奇点。包装函数为奇点返回 np.nan，以指示该点未定义被积数。如果需要，**kwargs 参数允许您将其他参数传递给四元函数。

通过使用此方法，集成应该继续进行，而不会引发集成警告。让我知道是否有效。