From my SymPy output I have the matrix shown below, which I must integrate in 2D. Currently I am doing it element-wise as shown below. This method works but it gets too slow (for both sympy.mpmath.quad and scipy.integrate.dblquad) for my real case (in which A and its functions are much bigger (see edit below):
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]
I was considering doing it in one shot like, but I'm not sure if this is the way to go:
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:
The real case has been made available in this link. Just unzip and run shadmehri_2012.py (is the author from were this example was taken from: Shadmehri et al. 2012).
I've started a bounty of 50 for the one who can do the following:
m=15 and n=15 in the code), I managed up to m=7 and n=7 in 32-bitThe current timing can be summarized below(measured with m=3 and n=3). From there it can be seen that the numerical integration is the bottleneck.
build trial functions = 0%
evaluating differential equations = 2%
lambdifying k1 = 22%
integrating k1 = 74%
lambdifying and integrating k2 = 2%
extracting eigenvalues = 0%
Related questions: about lambdify
I think you can avoid the lambdification time by switching to numerical evaluation at a different stage of the calculation.
Namely, your calculation seems to be diagonal in the sense that k1 and k2 are both of the form k = g^T X g where X is some 5x5 matrix (with differential ops inside, but that doesn't matter), and g is 5xM, with M large. Therefore k[i,j] = g.T[i,:] * X * g[:,j].
So you can just replace
for j in xrange(1,n+1):
for i in xrange(1,m+1):
g1 += [uu(i,j,x,t), 0, 0, 0, 0]
g2 += [ 0,vv(i,j,x,t), 0, 0, 0]
g3 += [ 0, 0,ww(i,j,x,t), 0, 0]
g4 += [ 0, 0, 0,bx(i,j,x,t), 0]
g5 += [ 0, 0, 0, 0,bt(i,j,x,t)]
g = Matrix( [g1, g2, g3, g4, g5] )
with
i1 = Symbol('i1')
j1 = Symbol('j1')
g1 = [uu(i1,j1,x,t), 0, 0, 0, 0]
g2 = [ 0,vv(i1,j1,x,t), 0, 0, 0]
g3 = [ 0, 0,ww(i1,j1,x,t), 0, 0]
g4 = [ 0, 0, 0,bx(i1,j1,x,t), 0]
g5 = [ 0, 0, 0, 0,bt(i1,j1,x,t)]
g_right = Matrix( [g1, g2, g3, g4, g5] )
i2 = Symbol('i2')
j2 = Symbol('j2')
g1 = [uu(i2,j2,x,t), 0, 0, 0, 0]
g2 = [ 0,vv(i2,j2,x,t), 0, 0, 0]
g3 = [ 0, 0,ww(i2,j2,x,t), 0, 0]
g4 = [ 0, 0, 0,bx(i2,j2,x,t), 0]
g5 = [ 0, 0, 0, 0,bt(i2,j2,x,t)]
g_left = Matrix( [g1, g2, g3, g4, g5] )
and
tmp = evaluateExpr( B*g ) k1 = r*tmp.transpose() * F * tmp k2 = r*g.transpose()*evaluateExpr(Bc*g) k2 = evaluateExpr( k2 )
by
tmp_right = evaluateExpr( B*g_right ) tmp_left = evaluateExpr( B*g_left ) k1 = r*tmp_left.transpose() * F * tmp_right k2 = r*g_left.transpose()*evaluateExpr(Bc*g_right) k2 = evaluateExpr( k2 )
Didn't test (past am), but you get the idea.
Now, instead of having a huge symbolic matrix which makes everything slow, you have two matrix indices for the trial function indices, and free parameters i1,j1 and i2,j2 which play their role and you should substitute integers into them in the end.
Since the matrix to lambdify is only 5x5, and needs to be lambdified only once outside all loops, the lambdification and simplification overhead is gone. Moreover, the problem fits easily into memory even for large m, n.
The integration is not any faster, but since the expressions are very small, you can easily e.g. dump them in Fortran or do something else smart.