scipy.quad trouble for decreasing functions over large ranges
I have a problem with scipy.quad . In short, I have a really long and complicated set of nested functions and integrals which include an integral of a decreasing function which must be integrated over the specific range 10^2 < x < 10^20.
To demonstrate this problem simply, consider the integral of y=x^(-2) between these values using numpy.quad :
import numpy as np
from scipy.integrate import quad
def func(x,z):
"decreasing function to test"
#print x,x**-2.
return x**(-2.)
#Test quad:
print 'Quad', quad(lambda x: func(x,0.),1e2,1e20)[0]
The true answer to this integral is 0.01, however quad returns a very small value, and by uncommenting the print line in the function you can see that it only integrates over the largest x values (or is biased that way due to the log scale of the x-axis). I need to find some way of fixing this problem.
I know that you can get the correct answer by using other methods such as the Simpsons or trapezoidal rules:
from scipy.integrate import simps
from numpy import trapz
xarray=np.logspace(2,20,1000)
print 'Simpson logspace x',simps(func(xarray,0.),xarray)
print 'Trapezoid logspace x',trapz(func(xarray,0.),xarray)
However these include passing an array into the function, which is not possible with my actual code. Thus I'd have to include a for loop to generate ay array to integrate with, which slows the whole program down unacceptably.
Is there any trick to making quad work for this sort of x range, or anything which works like quad which I can use instead?
链接地址: http://www.djcxy.com/p/69966.html上一篇: 在Python中读取特定的行序列