Reputation: 37
TypeError when I try to implement a function on a grid
Introduction
I have a pair of functions called numerical(x,y,z) and analytic(x,y,z), where numerical(x,y,z) is an approximation to analytic(x,y,z). I want to make some contour plots to see how these functions look in the x-y plane, so for the functions below I have set z=0. I have then tried to define the functions on a grid;
def numerical_plane(x,y):
return numerical(x, y, 0)
def analytic_plane(x,y):
R = np.sqrt(x**2+y**2)
return -G*M/(np.sqrt(R**2 + (a+np.sqrt(b**2))**2))
x = np.linspace(-20, 20, 50)
y = np.linspace(-20, 20, 50)
X, Y = np.meshgrid(x, y)
Z1 = numerical(X, Y)
Z2 = analytic(X, Y)
the function numerical(x,y,z) is defined at the bottom of this question.
Issue
Z2 works fine, but Z1 raises
TypeError: only size-1 arrays can be converted to Python scalars
As this is a TypeError, I tried checking the type returned by numerical(1,2); it was a float. analytic(1,2) on the other hand returns a numpy.float64. I tried forcing numerical(x,y) to return a numpy.float64 but this didn't actually help, so maybe I'm on the wrong track. Does anyone understand what's going wrong here?
If you need more information about the function numerical(x,y,z), then I've written it's definition here
def rho(x,y,z,M = 5e10,a = 4,b =0.8):
R = np.sqrt(x**2+y**2)
rho = ((b**2*M)/(4*np.pi))*(a*R**2+(a+3*np.sqrt(z**2+b**2))*(a+np.sqrt(z**2+b**2))**2)\
/((R**2+(a+np.sqrt(z**2+b**2))**2)**(2.5)*(z**2+b**2)**(1.5))
return rho
def numerical(x, y, z, s=0.01):
integrand = lambda xx, yy, zz: -G*rho(xx,yy,zz,M,a,b)/(np.sqrt((xx- x)**2+(yy-y)**2+(zz-z)**2+s**2))
phi = nquad(integrand, ranges = [(-np.inf, np.inf), (-np.inf, np.inf), (-np.inf, np.inf)], opts = {'epsrel' : 1e-1})[0]
return phi
It involves calling another function rho(x,y,z, ...), and then performing an integration over all space.
Here is also the full-traceback ...
TypeError Traceback (most recent call last)
<ipython-input-62-cb38a0d34db9> in <module>()
11 X, Y = np.meshgrid(x, y)
12 Z1 = analytic(X, Y)
---> 13 Z2 = numeric(X, Y)
<ipython-input-62-cb38a0d34db9> in numeric(x, y)
1 def numeric(x,y):
----> 2 return numerical_phi_MN(x, y, 0)
3
4 def analytic(x,y):
5 R = np.sqrt(x**2+y**2)
<ipython-input-43-c7e02143837e> in numerical_phi_MN(x, y, z, s)
14
15 integrand = lambda xx, yy, zz: -G*rho_MN(xx,yy,zz,M,a,b)/(np.sqrt((xx-x)**2+(yy-y)**2+(zz-z)**2+s**2))
---> 16 phi = nquad(integrand, ranges = [(-np.inf, np.inf), (-np.inf, np.inf), (-np.inf, np.inf)], opts = {'epsrel' : 1e-1})[0]
17
18 return phi
/opt/python-2.7.14/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc in nquad(func, ranges, args, opts, full_output)
712 else:
713 opts = [opt if callable(opt) else _OptFunc(opt) for opt in opts]
--> 714 return _NQuad(func, ranges, opts, full_output).integrate(*args)
715
716
/opt/python-2.7.14/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc in integrate(self, *args, **kwargs)
767 f = partial(self.integrate, depth=depth+1)
768 quad_r = quad(f, low, high, args=args, full_output=self.full_output,
--> 769 **opt)
770 value = quad_r[0]
771 abserr = quad_r[1]
/opt/python-2.7.14/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
321 if (weight is None):
322 retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
--> 323 points)
324 else:
325 retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel,
/opt/python-2.7.14/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
388 return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
389 else:
--> 390 return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
391 else:
392 if infbounds != 0:
/opt/python-2.7.14/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc in integrate(self, *args, **kwargs)
767 f = partial(self.integrate, depth=depth+1)
768 quad_r = quad(f, low, high, args=args, full_output=self.full_output,
--> 769 **opt)
770 value = quad_r[0]
771 abserr = quad_r[1]
/opt/python-2.7.14/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
321 if (weight is None):
322 retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
--> 323 points)
324 else:
325 retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel,
/opt/python-2.7.14/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
388 return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
389 else:
--> 390 return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
391 else:
392 if infbounds != 0:
/opt/python-2.7.14/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc in integrate(self, *args, **kwargs)
767 f = partial(self.integrate, depth=depth+1)
768 quad_r = quad(f, low, high, args=args, full_output=self.full_output,
--> 769 **opt)
770 value = quad_r[0]
771 abserr = quad_r[1]
/opt/python-2.7.14/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
321 if (weight is None):
322 retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
--> 323 points)
324 else:
325 retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel,
/opt/python-2.7.14/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
388 return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
389 else:
--> 390 return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
391 else:
392 if infbounds != 0:
TypeError: only size-1 arrays can be converted to Python scalars
Many thanks
Upvotes: 0
Views: 38
Answers (1)
Reputation: 7186
The problem is that scipy.integrate.nquad can only integrate a single (n-dimensional) function. Therefore it will always return a single value, instead of applying to a multidimensional vector element-wise.
You can get around this by making the function vectorized, using numpy.vectorize. However, note what is being said in the documentation:
The vectorize function is provided primarily for convenience, not for performance. The implementation is essentially a for loop.
With this (and G = 1, M = 1, rho_MN = rho), it runs (albeit taking a long time, will let you know when it finishes):
def numerical_plane(x, y):
return np.vectorize(numerical)(x, y, 0)
Upvotes: 1
Related Questions
- How to return a grid as a function?
- Passing numpy array to function giving type error
- What am I doing wrong in defining my function to get the error "numpy.ndarray" not callable?
- function gets function as an argument and throws 'numpy.ndarray' object is not callable
- Vecrtorized evluation of function defined by matrix over grid
- " numpy.ndarray is object is not callable " error in initiating a function
- TypeError: unsupported operand type(s) for -: 'float' and 'builtin_function_or_method' matrix
- TkInter grid function error message
- Python: how to evaluate a function on a grid
- Getting readable error on python while using grid function