title: scipy.integrate.quad on arrays of limits
body: <p>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. </p>
  
  <p>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.</p>
  
  <p>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?</p>
  
  <pre><code>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()
  """
  </code></pre>
  
  <p><strong>EDIT:</strong></p>
  
  <p>The actual function I need to integrate (<code>rho</code>) also contains two other factors. <code>rho0</code> is an array with the same length as <code>a_list</code> and <code>b_list</code>. <code>H</code> is a scalar.</p>
  
  <pre><code>def rho(x, rho0, H):
      return rho0 * np.exp(- x*x / (2*H*H))
  </code></pre>
  
  <p><strong>EDIT2:</strong></p>
  
  <p>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).</p>
  
  <pre><code>   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)
  </code></pre>