pythonnumpyscipycalculationnumerical-integration

Is it possible to somehow take this integral?


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))


Solution

  • 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))