Scipy分析性地积分分段函数

Question

在scipy中是否有针对分段函数的解析积分方法？例如，我有以下函数：

xrange_one, xrange_two = np.arange(0,4), np.arange(3,7)

part_one = lambda x: x + 3
part_two = lambda x: -2*x + 2

我想要集成这个分段函数的第一个矩：

func_one = lambda x: x * (x + 3)
func_two = lambda x: x * (-2*x + 2) 

是否有一种方法可以使用scipy integrate.quad或其他分析积分函数来执行此操作：

total = integrate.quad(func_one, 0, 3, func_two, 3, 6)

我不仅仅想把这两个部分分开集成。

- user2520932

---

“Analytical”和“scipy”不太兼容。如果您需要进行分析集成，请使用sympy。相反，scipy用于数值问题和数值积分。 - Andras Deak -- Слава Україні

“analytical”是什么意思？quad是数值积分，它不会先计算1/3*x**3 + 3/2*x**2表达式。np.polyint可以对这两个表达式执行不定多项式积分。 - hpaulj

好的，无论是分析法还是使用类似于quad这样的方法来传递要被积分的函数都可以。 - user2520932

Answer 1

Scipy不会为您执行解析积分，因为它是用于解决数值问题的。另一方面，Sympy可以精确处理简单的符号问题:

>>> import sympy as sym
>>> x = sym.symbols('x')
>>> f = sym.Piecewise((x*(x+3),x<3), (x*(-2*x+2),True))
>>> sym.integrate(f,(x,0,6))
-153/2

比较

>>> import scipy.integrate as integrate
>>> integrate.quad(lambda x:x*(x+3),0,3)[0] + integrate.quad(lambda x:x*(-2*x+2),3,6)[0]
-76.5
>>> -153/2.
-76.5

你可以先定义你的原始分段函数，然后将其乘以符号x ，然后对这个新函数进行解析积分。

另一个选择可能更接近你问题的精神，是数值定义分段函数，并在使用scipy之后。这仍然可以为您节省一些工作，但不会严格分析：

>>> f = lambda x: x*(x+3) if x<3 else x*(-2*x+2)
>>> integrate.quad(f,0,6)[0]
-76.5

这种方法中最完整的设置如下：

>>> f = lambda x: x+3 if x<3 else -2*x+2
>>> xf = lambda x: x*f(x)
>>> first_mom = integrate.quad(xf,0,6)[0]
>>> print(first_mom)
-76.5

首先，我们为f定义分段的lambda函数，然后将其与x相乘得到第一矩的被积函数。然后我们进行积分。

---

请注意，许多人不赞成将lambda绑定到变量上。如果您想正确地执行此操作，最好为您的分段函数定义一个命名函数，并仅在积分内部使用lambda（如果您不使用该被积函数，则可以忽略此建议）：

import scipy.integrate as integrate
def f(x):
    return x+3 if x<3 else -2*x+2

first_mom = integrate.quad(lambda x: x*f(x),0,6)[0]

- Andras Deak -- Слава Україні

Answer 2

在使用 numpy 的 poly 函数后，我得到了以下结果：

integrate.quad(lambda x:np.piecewise(x, [x < 3, x >= 3], 
     [lambda x: np.polyval([1,3,0],x), 
      lambda x: np.polyval([-2,2,0],x)]),
     0,6)

评估结果为：

(-76.5, 1.3489209749195652e-12)

有一个polyint函数来进行多项式积分。

In [1523]: np.polyint([1,3,0])
Out[1523]: array([ 0.33333333,  1.5       ,  0.        ,  0.        ])
In [1524]: np.polyint([-2,2,0])
Out[1524]: array([-0.66666667,  1.        ,  0.        ,  0.        ])

那是

x*(x+3) => x**2 + 3*x => np.poly1d([1,3,0]) => 1/3 x**3 + 3/2 x**2

所以，这两个 polyint 对象的适当端点差值是 analytical 解决方案：

In [1619]: np.diff(np.polyval(np.polyint([1,3,0]),[0,3])) +   
           np.diff(np.polyval(np.polyint([-2,2,0]),[3,6]))
Out[1619]: array([-76.5])

In [1621]: [np.polyval(np.polyint([1,3,0]),[0,3]), 
            np.polyval(np.polyint([-2,2,0]),[3,6])]
Out[1621]: [array([  0. ,  22.5]), array([  -9., -108.])]

- hpaulj

---

我甚至不知道np.piecewise是什么，太好了。 （我很高兴我不是唯一一个拥有考古学ipython历史记录的人。） - Andras Deak -- Слава Україні