utilizando SciPy para integrar una función que devuelve una matriz o matriz

Tengo una matriz simbólica que se puede express como:

from sympy import lambdify, Matrix g_sympy = Matrix([[ x, 2*x, 3*x, 4*x, 5*x, 6*x, 7*x, 8*x, 9*x, 10*x], [x**2, x**3, x**4, x**5, x**6, x**7, x**8, x**9, x**10, x**11]]) g = lambdify( (x), g_sympy ) 

De modo que para cada x obtengo una matriz diferente:

 g(1.) # matrix([[ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.], # [ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]]) g(2.) # matrix([[ 2.00e+00, 4.00e+00, 6.00e+00, 8.00e+00, 1.00e+01, 1.20e+01, 1.40e+01, 1.60e+01, 1.80e+01, 2.00e+01], # [ 4.00e+00, 8.00e+00, 1.60e+01, 3.20e+01, 6.40e+01, 1.28e+02, 2.56e+02, 5.12e+02, 1.02e+03, 2.05e+03]]) 

y así…

Necesito integrar numéricamente g sobre x , digamos from 0. to 100. (en el caso real, la integral no tiene una solución exacta) y en mi enfoque actual tengo que lambdify cada elemento en g e integrarlo individualmente. Estoy usando quad para hacer una integración de elementos como:

 ans = np.zeros( g_sympy.shape ) for (i,j), func_sympy in ndenumerate(g_sympy): func = lambdify( (x), func_sympy) ans[i,j] = quad( func, 0., 100. ) 

Hay dos problemas aquí: 1) lambdify usado muchas veces y 2) para loop ; y creo que el primero es el cuello de botella, porque la matriz g_sympy tiene un máximo de 10000 términos (lo que no es un gran problema para un bucle for).

Como se muestra arriba, lambdify permite la evaluación de toda la matriz, así que pensé: “¿Hay alguna forma de integrar toda la matriz?”

scipy.integrate.quadrature tiene un parámetro vec_func que me dio esperanza. Esperaba algo como:

 g_int = quadrature( g, x1, x2 ) 

para obtener la matriz totalmente integrada, pero proporciona el ValueError: matriz debe ser bidimensional


EDIT: Lo que estoy tratando de hacer se puede hacer aparentemente en Matlab usando quadv y ya ha sido discutido para SciPy

El caso real ha sido puesto a disposición aquí .

Para ejecutarlo necesitarás:

  • adormecido
  • scipy
  • matplotlib
  • simpy

Simplemente ejecute: "python curved_beam_mrs.py" .

Verá que el procedimiento ya es lento, principalmente debido a la integración, indicada por TODO en el archivo curved_beam.py .

Será mucho más lento si elimina el comentario indicado después del TODO en el archivo curved_beam_mrs.py .

La matriz de funciones que se integra se muestra en el archivo print.txt .

¡Gracias!

El primer argumento de quadrature o quadrature debe ser invocable. El argumento vec_func de la quadrature refiere a si el argumento de este llamable es un vector (posiblemente multidimensional). Técnicamente, puedes vectorize el propio quad :

 >>> from math import sin, cos, pi >>> from scipy.integrate import quad >>> from numpy import vectorize >>> a = [sin, cos] >>> vectorize(quad)(a, 0, pi) (array([ 2.00000000e+00, 4.92255263e-17]), array([ 2.22044605e-14, 2.21022394e-14])) 

Pero eso es simplemente equivalente a un bucle explícito sobre los elementos de a . Específicamente, no le dará ningún aumento de rendimiento, si eso es lo que está buscando. Entonces, en general, la pregunta es por qué y qué es exactamente lo que intentas lograr aquí.

Vectorización de las reglas integrales trapezoidales y simpson. Trapezoidal solo se copia y se pega de otro proyecto que usó espacio de registro en lugar de espacio lins para que pueda utilizar cuadrículas no uniformes.

 def trap(func,a,b,num): xlinear=np.linspace(a,b,num) slicevol=np.diff(xlinear) output=integrand(xlinear) output=output[:,:-1]+output[:,1:] return np.dot(output,slicevol)/2 def simpson(func,a,b,num): a=float(a) b=float(b) h=(ba)/num output=4*np.sum(integrand(a+h*np.arange(1,num,2)),axis=1) output+=2*np.sum(integrand(a+h*np.arange(2,num-1,2)),axis=1) output+=np.sum(integrand(b),axis=1) output+=np.sum(integrand(a),axis=1) return output*h/3 def integrand(rlin): first=np.arange(1,11)[:,None] second=np.arange(2,12)[:,None] return np.vstack((rlin*first,np.power(rlin,second))) 

Examine los errores relativos acumulativos trapazoidales y simpsons:

 b=float(100) first=np.arange(1,11)*(b**2)/2 second=np.power(b,np.arange(3,13))/np.arange(3,13) exact=np.vstack((first,second)) for x in range(3): num=x*100+100 tr=trap(integrand,0,b,num).reshape(2,-1) si=simpson(integrand,0,b,num).reshape(2,-1) rel_trap=np.sum(abs((tr-exact)/exact))*100 rel_simp=np.sum(abs((si-exact)/exact))*100 print 'Number of points:',num,'Trap Rel',round(rel_trap,6),'Simp Rel',round(rel_simp,6) Number of points: 100 Trap Rel 0.4846 Simp Rel 0.000171 Number of points: 200 Trap Rel 0.119944 Simp Rel 1.1e-05 Number of points: 300 Trap Rel 0.053131 Simp Rel 2e-06 

Cronométralo. Tenga en cuenta que ambas reglas trapezoidales utilizan 200 puntos, mientras que los simpsons están cronometrados en solo 100 según la convergencia anterior. Lo siento, no tengo sympy:

 s=""" import numpy as np from scipy.integrate import trapz def integrand(rlin): first=np.arange(1,11)[:,None] second=np.arange(2,12)[:,None] return np.vstack((rlin*first,np.power(rlin,second))) def trap(func,a,b,num): xlinear=np.linspace(a,b,num) slicevol=np.diff(xlinear) output=integrand(xlinear) output=output[:,:-1]+output[:,1:] return np.dot(output,slicevol)/2 def simpson(func,a,b,num): a=float(a) b=float(b) h=(ba)/num output=4*np.sum(integrand(a+h*np.arange(1,num,2)),axis=1) output+=2*np.sum(integrand(a+h*np.arange(2,num-1,2)),axis=1) output+=np.sum(integrand(b),axis=1) output+=np.sum(integrand(a),axis=1) return output*h/3 def simpson2(func,a,b,num): a=float(a) b=float(b) h=(ba)/num p1=a+h*np.arange(1,num,2) p2=a+h*np.arange(2,num-1,2) points=np.hstack((p1,p2,a,b)) mult=np.hstack((np.repeat(4,p1.shape[0]),np.repeat(2,p2.shape[0]),1,1)) return np.dot(integrand(points),mult)*h/3 def x2(x): return x**2 def x3(x): return x**3 def x4(x): return x**4 def x5(x): return x**5 def x5(x): return x**5 def x6(x): return x**6 def x7(x): return x**7 def x8(x): return x**8 def x9(x): return x**9 def x10(x): return x**10 def x11(x): return x**11 def xt5(x): return 5*x """ zhenya=""" a=[xt5,xt5,xt5,xt5,xt5,xt5,xt5,xt5,xt5,xt5,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11] vectorize(quad)(a, 0, 100) """ usethedeathstar=""" g=lambda x: np.array([[x,2*x,3*x,4*x,5*x,6*x,7*x,8*x,9*x,10*x],[x**2,x**3,x**4,x**5,x**6,x**7,x**8,x**9,x**10,x**11]]) xv=np.linspace(0,100,200) trapz(g(xv)) """ vectrap=""" trap(integrand,0,100,200) """ vecsimp=""" simpson(integrand,0,100,100) """ vecsimp2=""" simpson2(integrand,0,100,100) """ print 'zhenya took',timeit.timeit(zhenya,setup=s,number=100),'seconds.' print 'usethedeathstar took',timeit.timeit(usethedeathstar,setup=s,number=100),'seconds.' print 'vectrap took',timeit.timeit(vectrap,setup=s,number=100),'seconds.' print 'vecsimp took',timeit.timeit(vecsimp,setup=s,number=100),'seconds.' print 'vecsimp2 took',timeit.timeit(vecsimp2,setup=s,number=100),'seconds.' 

Resultados:

 zhenya took 0.0500509738922 seconds. usethedeathstar took 0.109386920929 seconds. vectrap took 0.041011095047 seconds. vecsimp took 0.0376999378204 seconds. vecsimp2 took 0.0311458110809 seconds. 

Algo para señalar en los tiempos es que la respuesta de zhenya debería ser mucho más precisa. Creo que todo es correcto, por favor avíseme si se requieren cambios.

Si proporciona las funciones y el rango que utilizará, probablemente pueda preparar algo mejor para su sistema. ¿También estaría interesado en utilizar núcleos / nodos adicionales?

En el caso real, la integral no tiene una solución exacta, ¿te refieres a las singularidades? ¿Podría ser más preciso en esto, así como en el tamaño de la matriz que desea integrar? Debo admitir que Sympy es terriblemente lento cuando se trata de algunas cosas (no estoy seguro de que la integración sea parte de ella, pero prefiero mantenerme alejado de Sympy y apegarme a la solución numpy). ¿Desea obtener una solución más elegante, al hacerlo con una matriz o una más rápida?

-nota: aparentemente no puedo agregar comentarios a tu publicación para pedir estas cosas, así que tuve que publicar esto como respuesta, tal vez esto se deba a que no tengo suficiente reputación, ¿o no? –

editar: algo como esto?

  import numpy from scipy.integrate import trapz g=lambda x: numpy.array([[x,2*x,3*x],[x**2,x**3,x**4]]) xv=numpy.linspace(0,100,200) print trapz(g(xv)) 

habiendo visto que quieres integrar cosas como sum (a * sen (bx + c) ^ n * cos (dx + e) ​​^ m), para diferentes coeficientes para a, b, c, d, e, m, n , sugiero hacer todo eso analíticamente. (debería tener alguna fórmula para eso, ya que simplemente puedes reescribir sin a exponenciales complejos

Otra cosa que noté cuando verifiqué esas funciones un poco mejor, es que sin (a * x + pi / 2) y sin (a * x + pi) y cosas por el estilo se pueden reescribir en cos o sin de una manera que elimine el pi / 2 o pi. También lo que veo es que solo mirando el primer elemento en su matriz de funciones:

 a*sin(bx+c)^2+d*cos(bx+c)^2 = a*(sin^2+cos^2)+(da)*cos(bx+c)^2 = a+(da)*cos(bx+c)^2 

lo que también simplifica los cálculos. Si tenía las fórmulas de una manera que no involucraba un archivo txt masivo, verifique cuál es la fórmula más general que necesita integrar, pero creo que es algo así como un * sin ^ n (bx + c) * cos ^ m (dx + e), donde myn son 0 1 o 2, y esas cosas se pueden simplificar en algo que se puede integrar analíticamente. Entonces, si descubre la función analítica más general que tiene, puede hacer algo como:

 f=lambda x: [[s1(x),s2(x)],[s3(x),s4(x)]] res=f(x2)-f(x1) 

donde s1 (x) etc son solo las versiones analíticamente integradas de sus funciones?

(Realmente no estoy planeando revisar todo el código para ver lo que hace el rest, pero ¿simplemente está integrando esas funciones en el archivo txt de la A a la B o algo así? ¿O hay algún lugar como ese que toma el cuadrado de ¿Cada función o cualquier cosa que pueda arruinar la posibilidad de hacerlo analíticamente?)

Esto debería simplificar sus integrales, supongo?

primera integral y: segunda

hmm, ese segundo enlace no funciona, pero entiendes la idea del primero, supongo

Edite, ya que no desea soluciones analíticas: la mejora sigue siendo deshacerse de sympy:

 from sympy import sin as SIN from numpy import sin as SIN2 from scipy.integrate import trapz import time import numpy as np def integrand(rlin): first=np.arange(1,11)[:,None] second=np.arange(2,12)[:,None] return np.vstack((rlin*first,np.power(rlin,second))) def simpson2(func,a,b,num): a=float(a) b=float(b) h=(ba)/num p1=a+h*np.arange(1,num,2) p2=a+h*np.arange(2,num-1,2) points=np.hstack((p1,p2,a,b)) mult=np.hstack((np.repeat(4,p1.shape[0]),np.repeat(2,p2.shape[0]),1,1)) return np.dot(integrand(points),mult)*h/3 A=np.linspace(0,100.,200) B=lambda x: SIN(x) C=lambda x: SIN2(x) t0=time.time() D=simpson2(B,0,100.,200) print time.time()-t0 t1=time.time() E=trapz(C(A)) print time.time()-t1 t2=time.time() F=simpson2(C,0,100.,200) print time.time()-t2 

resultados en:

 0.000764131546021 sec for the faster method, but when using sympy 7.58171081543e-05 sec for my slower method, but which uses numpy 0.000519037246704 sec for the faster method, when using numpy, 

conclusión: use numpy, ditch sympy, (mi método de numpy más lento es en realidad más rápido en este caso, porque en este ejemplo solo lo probé en una función sin, en lugar de en un ndarray de ellos, pero el punto de desechar simpy aún permanece cuando se compara la hora de la versión numpy del método más rápido con la de la versión sympy del método más rápido)

Podría haber encontrado alguna forma interesante de hacer esto, a costa de definir diferentes símbolos para la matriz g_symp :

 import numpy as np from scipy.integrate import quad import sympy as sy @np.vectorize def vec_lambdify(var, expr, *args, **kw): return sy.lambdify(var, expr, *args, **kw) @np.vectorize def vec_quad(f, a, b, *args, **kw): return quad(f, a, b, *args, **kw)[0] Y = sy.symbols("y1:11") x = sy.symbols("x") mul_x = [y.subs(y,x*(i+1)) for (i,y) in enumerate(Y)] pow_x = [y.subs(y,x**(i+1)) for (i,y) in enumerate(Y)] g_sympy = np.array(mul_x + pow_x).reshape((2,10)) X = x*np.ones_like(g_sympy) G = vec_lambdify(X, g_sympy) I = vec_quad(G, 0, 100) print(I) 

con resultados:

 [[ 5.00000000e+03 1.00000000e+04 1.50000000e+04 2.00000000e+04 2.50000000e+04 3.00000000e+04 3.50000000e+04 4.00000000e+04 4.50000000e+04 5.00000000e+04] [ 5.00000000e+03 3.33333333e+05 2.50000000e+07 2.00000000e+09 1.66666667e+11 1.42857143e+13 1.25000000e+15 1.11111111e+17 1.00000000e+19 9.09090909e+20]] 

y usando la magia ipython %timeit vec_quad(G,0,100) obtuve

 1000 loops, best of 3: 527 µs per loop 

Creo que este enfoque es algo más limpio, a pesar de los malabares con los símbolos.

Quadpy (un proyecto mío) hace cuadratura vectorizada. Esta

 import numpy import quadpy def f(x): return [ [x, 2*x, 3*x, 4*x, 5*x, 6*x, 7*x, 8*x, 9*x, 10*x], [x**2, x**3, x**4, x**5, x**6, x**7, x**8, x**9, x**10, x**11] ] sol = quadpy.line_segment.integrate( f, numpy.array([[0.0], [100.0]]), quadpy.line_segment.GaussLegendre(6) ) print(sol) 

da

 [[[ 5.00000000e+03] [ 1.00000000e+04] [ 1.50000000e+04] [ 2.00000000e+04] [ 2.50000000e+04] [ 3.00000000e+04] [ 3.50000000e+04] [ 4.00000000e+04] [ 4.50000000e+04] [ 5.00000000e+04]] [[ 3.33333333e+05] [ 2.50000000e+07] [ 2.00000000e+09] [ 1.66666667e+11] [ 1.42857143e+13] [ 1.25000000e+15] [ 1.11111111e+17] [ 1.00000000e+19] [ 9.09090909e+20] [ 8.33333333e+22]]] 

Tiempos:

 %timeit quadpy.line_segment.integrate(f, [0.0, 100.0], gl) 10000 loops, best of 3: 65.1 µs per loop