Integración numérica sobre una matriz de funciones, SymPy y SciPy

Desde mi salida de SymPy tengo la matriz que se muestra a continuación, que debo integrar en 2D. Actualmente lo estoy haciendo de forma elemental como se muestra a continuación. Este método funciona pero se vuelve muy lento (para ambossympy.mpmath.quad yscipy.integrate.dblquad) para mi caso real (en el queA y sus funciones son mucho más grandes (ver edición abajo):

from sympy import Matrix, sin, cos
import sympy
import scipy
sympy.var( 'x, t' )
A = Matrix([[(sin(2-0.1*x)*sin(t)*x+cos(2-0.1*x)*cos(t)*x)*cos(3-0.1*x)*cos(t)],
            [(cos(2-0.1*x)*sin(t)*x+sin(2-0.1*x)*cos(t)*x)*sin(3-0.1*x)*cos(t)],
            [(cos(2-0.1*x)*sin(t)*x+cos(2-0.1*x)*sin(t)*x)*sin(3-0.1*x)*sin(t)]])

# integration intervals
x1,x2,t1,t2 = (30, 75, 0, 2*scipy.pi)

# element-wise integration
from sympy.utilities import lambdify
from sympy.mpmath import quad
from scipy.integrate import dblquad
A_int1 = scipy.zeros( A.shape, dtype=float )
A_int2 = scipy.zeros( A.shape, dtype=float )
for (i,j), expr in scipy.ndenumerate(A):
    tmp = lambdify( (x,t), expr, 'math' )
    A_int1[i,j] = quad( tmp, (x1, x2), (t1, t2) )
    # or (in scipy)
    A_int2[i,j] = dblquad( tmp, t1, t2, lambda x:x1, lambda x:x2 )[0]

Estaba considerando hacerlo de una sola vez, pero no estoy seguro de que este sea el camino a seguir:

A_eval = lambdify( (x,t), A, 'math' )
A_int1 = sympy.quad( A_eval, (x1, x2), (t1, t2)                 
# or (in scipy)
A_int2 = scipy.integrate.dblquad( A_eval, t1, t2, lambda x: x1, lambda x: x2 )[0]

EDITAR: El caso real ha sido puesto a disposición eneste enlace. Solo descomprime y correshadmehri_2012.py (es el autor de donde se tomó este ejemplo de:Shadmehri et al. 2012). Comencé una recompensa de 50 para el que puede hacer lo siguiente:

hacerlo razonablemente más rápido que la pregunta propuesta
administrar para ejecutar sin dar error de memoria incluso con una serie de términosm=15 yn=15 en el código), logré hastam=7 yn=7 en 32 bits

El tiempo actual se puede resumir a continuación (medido con m = 3 y n = 3). Desde allí se puede ver que la integración numérica es el cuello de botella.

construir funciones de prueba = 0%
evaluando ecuaciones diferenciales = 2%
lambdifying k1 = 22%
integrando k1 = 74%
lambdificación e integración k2 = 2%
extracción de valores propios = 0%

Preguntas relacionadas:sobre lambdify

Respuestas a la pregunta(2)

Su respuesta a la pregunta