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