CODEWITHSUNDEEP
pythonarraysnumpyscipyintegral
PaulMag
PaulMag

Reputation: 4148

scipy.integrate.quad on arrays of limits

quad from scipy.integrate needs the arguments func, a, b. Where func is a function to integrate and a and b is the lower and upper integration limits, respectively. a and b has to be numbers.

I have a situation where I need to evaluate the integral of a function for hundreds of thousands of different a, b and sum the results. This takes a long time to loop through. I tried to just give quad arrays for a and b, hoping quad would return the corresponding array, but that didn't work.

Here is a code illustrating what I'm trying to do, with both the Python loop that works, but is very slow, and my attempt at vectorizing that doesn't work. Any suggestions on how to this problem can be solved in a fast (numpy-is) way?

import numpy as np
from scipy.integrate import quad

# The function I need to integrate:
def f(x):
    return np.exp(-x*x)

# The large lists of different limits:
a_list = np.linspace(1, 100, 1e5)
b_list = np.linspace(2, 200, 1e5)

# Slow loop:
total_slow = 0  # Sum of all the integrals.
for a, b in zip(a_list, b_list):
    total_slow += quad(f, a, b)[0]
        # (quad returns a tuple where the first index is the result,
        # therefore the [0])

# Vectorized approach (which doesn't work):
total_fast = np.sum(quad(f, a_list, b_list)[0])

"""This error is returned:

 line 329, in _quad
    if (b != Inf and a != -Inf):
ValueError: The truth value of an array with more than 
one element is ambiguous. Use a.any() or a.all()
"""

EDIT:

The actual function I need to integrate (rho) also contains two other factors. rho0 is an array with the same length as a_list and b_list. H is a scalar.

def rho(x, rho0, H):
    return rho0 * np.exp(- x*x / (2*H*H))

EDIT2:

Profiling of different solutions. ´space_sylinder` is the function where the integral happens. Warren Weckesser suggestion is as fast as passing the arrays through a simple analytical function and ~500 times faster than the slow Python loop (note by the number of calls that the program didn't even finish and it still used 657 seconds).

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
# Analytic (but wrong) approximation to the solution of the integral:
      108    1.850    0.017    2.583    0.024 DensityMap.py:473(space_sylinder)
# Slow python loop using scipy.integrate.quad:
       69   19.223    0.279  657.647    9.531 DensityMap.py:474(space_sylinder)
# Vectorized scipy.special.erf (Warren Weckesser's suggestion):
      108    1.786    0.017    2.517    0.023 DensityMap.py:475(space_sylinder)

Upvotes: 2

Views: 3461

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114781

The integral of exp(-x*x) is a scaled version of the error function, so you can use scipy.special.erf to compute the integral. Given scalars a and b, the integral of your function from a to b is 0.5*np.sqrt(np.pi)*(erf(b) - erf(a)).

erf is a "ufunc", which means it handles array arguments. Given a_list and b_list, your calculation can be written as

total = 0.5*np.sqrt(np.pi)*(erf(b_list) - erf(a_list)).sum()

The function rho can also be handled with erf, by using the appropriate scaling:

g = np.sqrt(2)*H
total = g*rho0*0.5*np.sqrt(np.pi)*(erf(b_list/g) - erf(a_list/g)).sum()

Check this against your slow solution before relying on it. For some values, the subtraction of the erf functions can result in a significant loss of precision.

Upvotes: 3

Related Questions