I have the following code, but it shows error:
IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties. If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges. Perhaps a special-purpose integrator should be used.
potC=sc.integrate.quad(lambda r: Psi(r,n2)*(-1/r)Psi(r,n1)(r**2),0,np.inf)
How to fix it?
import scipy as sc
import numpy as np
def Psi(r,n):
return 2*np.exp(-r/n)*np.sqrt(n)*sc.special.hyp1f1(1-n, 2, 2*r/n)/n**2
def p(n1,n2):
potC=sc.integrate.quad(lambda r: Psi(r,n2)*(-1/r)*Psi(r,n1)*(r**2),0,np.inf)
pot=potC[0]
return pot
print(p(15,15))
The error literally says what your problem is. Your function is not "well behaved" in some regions.
For example with n = 15 and r = 50 your special.hyp1f1(-14, 2, 2*50/15) result in NaN. I am not familiar with this function, so I do not know if this is the expected behaviour, but this is what happens.
You can try to isolate these points and exclude them from the integration (if you know lower and upper bounds of the function's value (if it is defined) you can also update the expected error of your integration) and just integrate in the well behaved regions.
If it is a bug in scipy, then please report it to them.
Ps.: Tested with scipy 1.8.0
Ps2.: With some reading I found, that you can get the values correctly, if you do your calculations with complex number, so the following code gives you a value:
import scipy as sc
from scipy import integrate
from scipy import special
import numpy as np
def Psi(r,n):
r = np.array(r,dtype=complex)
return 2*np.exp(-r/n)*np.sqrt(n)*special.hyp1f1(1-n, 2, 2*r/n)/n**2
def p(n1,n2):
potC=integrate.quad(lambda r: Psi(r,n2)*(-1/r)*Psi(r,n1)*(r**2),0,np.inf)
pot=potC[0]
return pot
print(p(15,15))