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