Integration with "quad" and "quadrature" from Python/SciPy
After reading this and that, it occurs to me that both "quad" and "quadrature" should be interchangeable*, atleast syntax-wise. Strangely it does seem they are not:
from scipy.integrate import quad as q
#from scipy.integrate import quadrature as q
def myfunc(x):
return x
def integr():
return q(myfunc, 0, 1)[0]
print integr()
def myfunc2(x, y):
return x + y
def integr2(y):
return q(myfunc2, 0, 1, args=(y))[0]
#return q(myfunc2, 0, 1, args=[y])[0]
print integr2(10)
... the example runs fine for "quad", but not for "quadrature" - I end up with:
Traceback (most recent call last):
File "./test.py", line 38, in <module>
print integr2(10)
File "./test.py", line 36, in integr2
return q(myfunc2, 0, 1, args=(y))[0]
File "/usr/lib/python2.6/dist-packages/scipy/integrate/quadrature.py", line 136, in quadrature
newval = fixed_quad(vfunc, a, b, (), n)[0]
File "/usr/lib/python2.6/dist-packages/scipy/integrate/quadrature.py", line 48, in fixed_quad
return (b-a)/2.0*sum(w*func(y,*args),0), None
File "/usr/lib/python2.6/dist-packages/scipy/integrate/quadrature.py", line 77, in vfunc
return func(x, *args)
TypeError: myfunc2() argument after * must be a sequence, not int
I have to switch the args tuple to a list (cf. commented line in integr2) even though the documentation says it should be a tuple. It seemed this is what the interpreter complains about ... (right?)
Is this intended? Or am I doing something wrong? In the end I'd like to be able to choose integration methods afterwards without having to change too much of the rest of the code.
*Actually开发者_如何转开发 I don't really get how to choose between the two. I do understand the difference between Gaussian quadrature and adaptive quadrature, but I don't know what "adaptive Gaussian quadrature" is supposed to mean - is the number of nodes adapted, if so how!?
The problem is in the line return q(myfunc2, 0, 1, args=(y))[0]
, specifically in the args=(y)
part. What you want is args=(y,)
(notice the comma after y
) or args=[y]
.
The issue is that in Python tuples are created with commas, not with parentheses. Look:
>>> a = (1,)
>>> b = (1)
>>> print a, type(a)
(1,) <type 'tuple'>
>>> print b, type(b)
1 <type 'int'>
精彩评论