Scipy analytically integrate piecewise function

Question

Is there a method in scipy for analytic integration of a piecewise function? For example I have:

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

I'd like to integrate the first moment of this piecewise function:

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

Is there a way with scipy integrate.quad or some other analytical integration function that do something like this:

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

I don't want to just integrate the two pieces separately.

Answer 1

After playing with the numpy poly functions, I came up with:

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)

evaluates to:

(-76.5, 1.3489209749195652e-12)

There is a polyint to do a polynomial integration

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

That is

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

So the analytical solution is the appropriate end point differences for these two polyint objects:

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

Answer 2

Scipy will not perform analytical integration for you, since it is made for solving numerical problems. Sympy, on the other hand, can handle simple symbolic problems exactly:

>>> 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

Compare

>>> 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

You could also define your original piecewise function first, then multiply it with the symbolic x, then integrate this new function analytically.

Another alternative, perhaps closer to the spirit of your question, might be to define the piecewise function numerically, and using scipy after all. This will still save you some work, but won't be strictly analytical:

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

The most complete setup with this approach:

>>> 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

First we define the piecewise lambda for f, then the integrand of the first moment, multiplying it with x. Then we do the integration.

---

Note that it is frowned upon by many to bind lambdas to variables. If you want to do this properly, you should probably define a named function for your piecewise function, and only use a lambda inside the integration (if otherwise you wouldn't use that integrand):

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]