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) 

我不想分开整合两件。

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，那么积第一时刻，与x相乘。然后我们进行整合。

---

请注意，它被许多人视为将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] 

Answer 2

摆弄后在numpypoly功能，我想出了：

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 

所以analytical解决方案是针对这两个polyint对象适当终点的差异：

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.])]