title: How to Lambdify a multivariate sympy expression for numeric integration efficiently
body: <p>I have this very complex symbolic function of three variables, f(x,y,z), which I need to evaluate over a point in z and integrate over all the domain in x and y.</p>
  <p>Summarizing, I've been using:</p>
  <pre><code>f_substituted=f.subs(z,valueforz)
  
  f_numeric=lambdify([x,y].f_substituted)
  
  #and then, integrate it a la Riemann (dblquad does not converge), by doing:
  
  N=1000
  q = np.linspace(-10, 10, N)
  p = np.linspace(-10, 10, N)
  dS=(q[1]-q[0])*(p[1]-p[0])
  q, p = np.meshgrid(q, p) #I meshgrid to have a 2D grid for integration.
  
  #I define a dataframe here because I have seen it is quicker than vectorizing the expression
  df = pd.DataFrame({'Q':np.concatenate(q),'P':np.concatenate(p)})
  integral=sum((rho_dot_00_x(df['Q'],df['P'])).dropna()).real*dS
  </code></pre>
  <p>My problem is that it works really slowly and I cannot make work ufuncify, or theano.</p>
  <p>Also, I tried to lambdify using &quot;numpy&quot; argument, but it ran way slowlier (I thought it should be quicker than math).
  I also tried evalf after the substitution, like this:</p>
  <pre><code>f_substituted=f.subs(z,valueforz).evalf()
  </code></pre>
  <p>but it was even slower.</p>
  <p>Any help to make the numeric evaluation quicker would be welcome.</p>
  <p>Thanks in advance.</p>