Integração Numérica sobre uma Matriz de Funções, SymPy e SciPy

Da minha saída do SymPy, tenho a matriz mostrada abaixo, que devo integrar em 2D. Atualmente estou fazendo o elemento como mostrado abaixo. Este método funciona, mas fica muito lento (para ambossympy.mpmath.quad escipy.integrate.dblquad) para o meu caso real (em queA e suas funções são muito maiores (veja abaixo):

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]

Eu estava pensando em fazer isso de uma vez, mas não tenho certeza se esse é o caminho 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]

EDIT: O caso real foi disponibilizado emesse link. Apenas descompacte e executeshadmehri_2012.py (é o autor de onde este exemplo foi retirado de:Shadmehri et al. 2012). Eu comecei uma recompensa de 50 para quem pode fazer o seguinte:

torná-lo razoavelmente mais rápido do que a questão proposta
consegue executar sem dar erro de memória, mesmo com um número de termosm=15 en=15 no código), consegui atém=7 en=7 em 32 bits

O tempo atual pode ser resumido abaixo (medido com m = 3 en = 3). A partir daí, pode-se ver que a integração numérica é o gargalo.

construir funções de teste = 0%
avaliando equações diferenciais = 2%
lambidificante k1 = 22%
integrando k1 = 74%
lambdifying e integrando k2 = 2%
extração de autovalores = 0%

Perguntas relacionadas:sobre lambdify

questionAnswers(2)

yourAnswerToTheQuestion