title: SciPy Can't Convert Complex to Float
body: <p>I am stuck at a problem, which I have simplified below. </p>
  
  <p>So, I have a function (matrix) which returns eigenvalues of a matrix.</p>
  
  <p>Second function (deriv) finds the derivatives of eigenvalues with respect to delta. SciPy's own derivative function is very slow, so I have used complex difference method.</p>
  
  <p>Then I take the double integral of integrand over r and theta. Instead of taking the double integral, I solve a differential equation of r, after taking one dimensional integral over theta, which makes the calculations faster. </p>
  
  <p>Finally, integral works. However, if I try to find the specific delta which satisfies the equation, it gives an error: Can't convert complex to float, which I don't understand where the imaginary number comes from. </p>
  
  <p>I am new to Python world, any help would be quite appreciated. Thank you.</p>
  
  <pre><code>from scipy.integrate import quad
  from numpy import linalg as LA
  import numpy as np
  from pylab import *
  
  from scipy.integrate import odeint
  from scipy.integrate import ode
  from scipy.optimize import fsolve
  
  
  #This routine is for the matrix and returns the eigenvalues
  def matrix(r, theta, delta1, delta2):
  
      mat = np.array([[r**2-1,r*np.cos(theta),0,delta1],r*np.cos(theta),r**2 - 1, -delta1,0],[0,-delta2,-r**2+1,-r*np.cos(theta)],[delta2,0,-r*np.cos(theta),-r**2+1]])
      return np.sort(LA.eigvals(mat))[2:4]
  
  #This routine takes the derivatives of eigenvalues with respect to the parameter delta. I set delta1 = delta2.
  def deriv(r, theta, delta):
  
      h = 0.00000000001
  
      return np.imag(matrix(r, theta, delta, delta+h*1j))/h
  
  #This is the integrand that we need to integrate over r and theta
  def integrand(theta,r, beta, delta):
      ee = matrix(r, theta, delta, delta)
      dd = deriv(r, theta, delta)
      return (np.tanh(0.5*beta*ee[0])*dd[0]+np.tanh(0.5*beta*ee[1])*dd[1])*r*r*np.sin(theta)
  
  #Just integrate over theta
  def polarint(y,r,beta,delta):
      return quad(integrand,0,np.pi,args = (r,beta,delta))[0]
  
  #Instead of integrating directly over r, solve the differential equation.
  #Lambda is the range of r.
  def givesclen(delta, beta, lam):
      y0 = 0
      t_out = np.linspace(0, lam, 2500);
      rt = odeint(polarint,y0,t_out,args=(beta,delta))
      temp = (rt[-1]/delta)/(4*np.pi**2)-t_out[-1]/(2*np.pi**2)
      return temp[0]
  
  #The goal is to find the delta; given sl, lam, beta
  #Such that the result of the integration is equal to sl
  def equationgap(delta, beta, lam,sl):
      return givesclen(delta, beta, lam)*4*np.pi - sl
  
  #Test if the equationgap works, the result should be close to zero!
  print equationgap(.5,40,500,.1744)
  
  #Now use the fsolve function should find the delta to be .5!
  #beta = 40
  #lam = 500
  #sl = 0.174
  #fsolve(equationgap,.6,args = (beta, lam, sl))
  </code></pre>
  
  <p>Edit: </p>
  
  <p>The error message is:</p>
  
  <pre><code>Traceback (most recent call last):
  File "test.py", line 38, in polarint
  return quad(integrand,0,np.pi,args = (r,beta,delta))[0]
  File "/usr/local/anaconda/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 281, in quad
  retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
  File "/usr/local/anaconda/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 345, in _quad
  return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
  File "test.py", line 30, in integrand
  dd = deriv(r, theta, delta)
  File "test.py", line 22, in deriv
  return np.imag(matrix(r, theta, delta, delta+h*1j))/h
  File "test.py", line 14, in matrix
  mat = np.array([[r**2-1,r*np.cos(theta),0,delta1],[r*np.cos(theta),r**2 - 1, -delta1,0],[0,-delta2,-r**2+1,-r*np.cos(theta)],[delta2,0,-r*np.cos(theta),-r**2+1]])
  TypeError: can't convert complex to float
  </code></pre>