I want to evaluate an integral but python is complaining that the integrand becomes large at some point during the integration.
The goal is to find the value of a variable called Sigma.
Sigma is the outcome of the integral, but Sigma also appears in the integrand of the integral and the problem is therefore solved by an iteration where you do an initial guess, then put the guess into the integral, evaluate the integral and then but the outcome back into the integral etc.
a, U0, E_F and t are real constants. Sigma can in general be an imaginary number. I then define the integrand with a function and do the iteration 100 times.
My code is as follows:
a = 1
U0 = 0.001
E_F = 0.10
t = taa/(a**2)
Sigma = 0
def integrand(k_x,a, U0, E_F, t,Sigma):
return a*U0**2/(24*np.pi) * 1/(E_F - 2*t*(1-np.cos(a*k_x))-Sigma)
Sigma = quad(integrand,-np.pi/a,np.pi/a, args=(a,U0,E_F,t,Sigma))
for i in range(0,100):
Sigma = quad(integrand,-np.pi/a,np.pi/a, args=(a,U0,E_F,t,Sigma[0]))
(taa is specified somewhere else but it's just a real constant).
When running the code, Python gives warnings: IntegrationWarning: The integral is probably divergent, or slowly convergent. Python does give an answer, but with an error about twice the value of the answer.
I'm supposed to add a small imaginary part to the integrand epsilon*i, which should make sure that the integral converges.
When you have a real integral, doesn't it matter that the you have a complex integrand?
I tried adjusting the code to:
a = 1
U0 = 0.001
E_F = 0.10
t = taa/(a**2)
Sigma = 0
epsilon = 0.01*1j
def integrand(k_x,a, U0, E_F, t,epsilon,Sigma):
return a*U0**2/(24*np.pi) * 1/(E_F +epsilon - 2*t*(1-np.cos(a*k_x))-Sigma)
Sigma = quad(integrand,-np.pi/a,np.pi/a, args=(a,U0,E_F,t,epsilon,Sigma))
for i in range(0,100):
Sigma = quad(integrand,-np.pi/a,np.pi/a, args=(a,U0,E_F,t,epsilon,Sigma[0]))
This way the integral should converge I think, but Python gives the error: Supplied function does not return a valid float.
Does someone know what I can do here? Should I use a different integration function? Should this integral theoretically converge this way?
(If someone is wondering about the nature of the problem: Sigma is the self-energy, the integral is over the Brillouin zone, E_F is the Fermi level and the small imaginary part has to do with Green's functions)
