scipy.integrate.quad en matrices de límites

quad from scipy.integrate necesita los argumentos func, a, b. Donde func es una función para integrar y a y b son los límites de integración inferior y superior, respectivamente. A y B tienen que ser números.

Tengo una situación en la que necesito evaluar la integral de una función para cientos de miles de a, b diferentes y sumr los resultados. Esto lleva mucho tiempo para pasar. Intenté simplemente dar matrices de quad para ayb, esperando que quad devolviera la matriz correspondiente, pero eso no funcionó.

Aquí hay un código que ilustra lo que estoy tratando de hacer, con el bucle de Python que funciona, pero es muy lento, y mi bash de vectorización no funciona. ¿Alguna sugerencia sobre cómo resolver este problema se puede resolver de manera rápida (sin tener en cuenta)?

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() """ 

EDITAR:

La función real que necesito integrar ( rho ) también contiene otros dos factores. rho0 es una matriz con la misma longitud que a_list y b_list . H es un escalar.

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

EDIT2:

Perfilado de diferentes soluciones. ´space_sylinder` es la función donde ocurre la integral. La sugerencia de Warren Weckesser es tan rápida como pasar las matrices a través de una función analítica simple y ~ 500 veces más rápida que el lento bucle de Python (observe por la cantidad de llamadas que el progtwig ni siquiera terminó y todavía se usó 657 segundos).

  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) 

La integral de exp(-x*x) es una versión escalada de la función de error , por lo que puede usar scipy.special.erf para calcular la integral. Dados los escalares a y b , la integral de su función de a a b es 0.5*np.sqrt(np.pi)*(erf(b) - erf(a)) .

erf es un “ufunc” , lo que significa que maneja argumentos de matriz. Dado a_list y b_list , su cálculo puede escribirse como

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

La función rho también se puede manejar con erf , usando la escala apropiada:

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

Compruebe esto contra su solución lenta antes de confiar en ella. Para algunos valores, la resta de las funciones erf puede resultar en una pérdida significativa de precisión.