scipynumerical-integrationintegralmpmathquad

"TypeError" message while performing quad integration method of scipy.integrate with mpmath functions


I am trying to calculate two integrals using scipy.integrate.quad. However, since gamma functions with negative first parameter are not defined in scipy, I had to choose the version from mpmath. After running the following code,

from scipy.integrate import *
from mpmath import *



low, up  = 5.630e5, 1.167e12
alpha, threshold = 1.05   , 2.15e10 
beta = 274

def g(x, beta, low, up):
    return gamma(-2/3) * (gammainc(-2/3, beta*(x/low)**3) - gammainc(-2/3, beta*(x/up)**3))

def Integrand1(x, low, threshold, alpha):
    return pow(x/threshold, alpha) * g

def Integrand2(x, up, threshold):
    return g

Integral1 = quad(Integrand1, low, threshold, args=(low, up, threshold, alpha, beta))
Integral2 = quad(Integrand2, threshold, up, args=(low, up, threshold, beta))

print(Integral1)
print(Integral2)

Here is the error message which I don't know how to handle and need help with it:

Traceback (most recent call last): File "test.py", line 19, in Integral1 = quad(Integrand1, low, threshold, args=(low, up, threshold, alpha, beta)) File "/home/username/anaconda3/lib/python3.6/site-packages/mpmath/calculus/quadrature.py", line 748, in quad points[0], prec, epsilon, m, verbose) File "/home/username/anaconda3/lib/python3.6/site-packages/mpmath/calculus/quadrature.py", line 215, in summation for i in xrange(len(points)-1): TypeError: object of type 'float' has no len()

I can only guess that the reason could be quad function is not compatible with integrals defined using mpmath.


Solution

  • Import statements

    Don't import * from two places, it's a recipe for name collision. MpMath has its own quad method, which replaces SciPy's quad in your code.

    from scipy.integrate import quad
    from mpmath import gamma, gammainc 
    

    Function parameters

    If you are calling function g, you have to provide arguments to it. So, write * g(x, beta, low, up) instead of * g.

    Of course, these arguments must also be available to the function that is calling g. Like this:

    def Integrand1(x, low, up, threshold, alpha, beta):
        return pow(x/threshold, alpha) * g(x, beta, low, up)
    
    def Integrand2(x, low, up, threshold, alpha, beta):
        return g(x, beta, low, up)
    
    Integral1 = quad(Integrand1, low, threshold, args=(low, up, threshold, alpha, beta))
    Integral2 = quad(Integrand2, threshold, up, args=(low, up, threshold, alpha, beta))
    

    Note that the arguments passed to Integrand functions match what they expect to receive. They get x, and everything listed in args parameter of quad.

    No errors are thrown by the code above. I am not sure the operation is meaningful mathematically, as you use threshold both for scaling and as an upper limit, but that's another story.