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.