Funciones de integración que devuelven una matriz en Python

Tengo una gran cantidad de datos para integrar y me gustaría encontrar una forma de hacerlo todo con solo matrices, y estaría dispuesto a comprometer la precisión para un aumento del rendimiento. Lo que tengo en mente es algo como esto:

import numpy import scipy a = np.array([1,2,3]) def func(x): return x**2 + x def func2(x): global a return a*x def integrand(x): return func(x)*func2(x) integrated = quad(integrand, 0, 1) 

Así que estoy tratando de integrar cada elemento en la matriz que sale de integrand .

Soy consciente de que existe la posibilidad de usar numpy.vectorize() siguiente manera:

 integrated = numpy.vectorize(scipy.integrate.quad)(integrand, 0, 1) 

Pero no puedo hacer que funcione. ¿Hay una manera de hacer esto en python?

Solución

Bueno, ahora que aprendí un poco más de python, puedo responder a esta pregunta si a alguien le resulta estable y tiene la misma pregunta. La forma de hacerlo es escribir las funciones como si fueran a tomar valores escalares, y no vectores como entradas. Así que siga mi código anterior, lo que tendríamos es algo como

 import numpy as np import scipy.integrate.quad a = np.array([1, 2, 3]) # arbitrary array, can be any size def func(x): return x**2 + x def func2(x, a): return a*x def integrand(x, a): return func(x)*func2(x, a) def integrated(a): integrated, tmp = scipy.integrate.quad(integrand, 0, 1, args = (a)) return integrated def vectorizeInt(): global a integrateArray = [] for i in range(len(a)): integrate = integrated(a[i]) integrateArray.append(integrate) return integrateArray 

No es que la variable que está integrando sobre debe ser la primera entrada a la función. Esto es necesario para scipy.integrate.quad. Si se está integrando a través de un método, es el segundo argumento después del self típico (es decir, x está integrado en def integrand(self, x, a): . También es necesario args = (a) para indicar a quad el valor de a en la función integrand . Si integrand tiene muchos argumentos, digamos def integrand(x, a, b, c, d): simplemente pones los argumentos en orden en args . Entonces eso sería args = (a, b, c, d) .

vectorize no proporcionará ayuda para mejorar el rendimiento del código que usa quad . Para usar quad , deberá llamarlo por separado para cada componente del valor devuelto por integrate .

Para una aproximación vectorizada pero menos precisa, puede usar numpy.trapz o scipy.integrate.simps .

La definición de su función (al menos la que se muestra en la pregunta) se implementa utilizando funciones numpy que admiten la difusión, por lo que, dada una cuadrícula de valores x en [0, 1], puede hacer esto:

 In [270]: x = np.linspace(0.0, 1.0, 9).reshape(-1,1) In [271]: x Out[271]: array([[ 0. ], [ 0.125], [ 0.25 ], [ 0.375], [ 0.5 ], [ 0.625], [ 0.75 ], [ 0.875], [ 1. ]]) In [272]: integrand(x) Out[272]: array([[ 0. , 0. , 0. ], [ 0.01757812, 0.03515625, 0.05273438], [ 0.078125 , 0.15625 , 0.234375 ], [ 0.19335938, 0.38671875, 0.58007812], [ 0.375 , 0.75 , 1.125 ], [ 0.63476562, 1.26953125, 1.90429688], [ 0.984375 , 1.96875 , 2.953125 ], [ 1.43554688, 2.87109375, 4.30664062], [ 2. , 4. , 6. ]]) 

Es decir, al hacer x una matriz con forma (n, 1), el valor devuelto por integrand(x) tiene forma (n, 3) . Hay una columna para cada valor en a .

Puede pasar ese valor a numpy.trapz() o scipy.integrate.simps() , utilizando axis=0 , para obtener las tres aproximaciones de las integrales. Probablemente querrás una cuadrícula más fina:

 In [292]: x = np.linspace(0.0, 1.0, 101).reshape(-1,1) In [293]: np.trapz(integrand(x), x, axis=0) Out[293]: array([ 0.583375, 1.16675 , 1.750125]) In [294]: simps(integrand(x), x, axis=0) Out[294]: array([ 0.58333333, 1.16666667, 1.75 ]) 

Compare eso con las llamadas repetidas a quad :

 In [296]: np.array([quad(lambda t: integrand(t)[k], 0, 1)[0] for k in range(len(a))]) Out[296]: array([ 0.58333333, 1.16666667, 1.75 ]) 

Su función de integrate (que asumo es solo un ejemplo) es un polinomio cúbico, para el cual la regla de Simpson da el resultado exacto. En general, no espere que los simps den una respuesta tan precisa.

Quadpy (un proyecto mío) está completamente vectorizado. Instalar con

 pip install quadpy 

y luego hacer

 import numpy import quadpy def integrand(x): return [numpy.sin(x), numpy.exp(x)] # ,... res = quadpy.line_segment.integrate( integrand, [0, 1], quadpy.line_segment.GaussLegendre(5) ) print(res) 

Salida:

 [ 0.45969769 1.71828183]