完整的追溯（截取的部分）表明该错误出在univariatespline对象的__call__方法中。因此，实际上的问题是，mpmath集成例程以其mpf小数形式输入，而scipy无法处理它们。

然后最简单的解决方法是将integrand参数的有问题的部分手动转换为浮点数：

integrand = lambda U: dxdu_u(float(U)) * mp.besselj(n,U)

通常，这容易产生数字错误（mpmath故意使用其高精度变量！），因此请谨慎操作。在这种特定情况下，这可能没问题，因为插值实际上是以双精度完成的。不过，最好检查一下结果。

一种可能的选择是避免使用mpmath并在weights中使用scipy.integrate.quad参数，请参见docs（向下滚动到weights="sin"部分）

另一种选择是始终坚持使用mpmath并在纯python中实现插值（这样，mpf对象可能很好，因为它们应该支持常用的算术）。一个简单的线性插值可能就足够了。如果不是，那么编写自己的三次样条插值器就没什么大不了了。

完整的追溯：

In [443]: f(0)                                                                  
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-443-6bfbdbfff9c4> in <module>
----> 1 f(0)

<ipython-input-440-7ebeff3611f6> in f(n)
      2     integrand = lambda U: dxdu_u(U) * mp.besselj(n,U)
      3     bjz = lambda nth: mp.besseljzero(n,nth)
----> 4     return mp.quadosc(integrand,[0,mp.inf],zeros=bjz)
      5 

/usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in quadosc(ctx,f,interval,omega,period,zeros)
    998         #    raise ValueError("zeros do not appear to be correctly indexed")
    999         n = 1
-> 1000         s = ctx.quadgl(f,[a,zeros(n)])
   1001         def term(k):
   1002             return ctx.quadgl(f,[zeros(k),zeros(k+1)])

/usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in quadgl(ctx,*args,**kwargs)
    807         """
    808         kwargs['method'] = 'gauss-legendre'
--> 809         return ctx.quad(*args,**kwargs)
    810 
    811     def quadosc(ctx,omega=None,period=None,zeros=None):

/usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in quad(ctx,*points,**kwargs)
    740             ctx.prec += 20
    741             if dim == 1:
--> 742                 v,err = rule.summation(f,points[0],prec,epsilon,m,verbose)
    743             elif dim == 2:
    744                 v,err = rule.summation(lambda x: \

/usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in summation(self,points,max_degree,verbose)
    230                     print("Integrating from %s to %s (degree %s of %s)" % \
    231                         (ctx.nstr(a),ctx.nstr(b),degree,max_degree))
--> 232                 results.append(self.sum_next(f,nodes,results,verbose))
    233                 if degree > 1:
    234                     err = self.estimate_error(results,epsilon)

/usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in sum_next(self,previous,verbose)
    252         case the quadrature rule is able to reuse them.
    253         """
--> 254         return self.ctx.fdot((w,f(x)) for (x,w) in nodes)
    255 
    256 

/usr/local/lib/python3.6/dist-packages/mpmath/ctx_mp_python.py in fdot(ctx,A,B,conjugate)
    942         hasattr_ = hasattr
    943         types = (ctx.mpf,ctx.mpc)
--> 944         for a,b in A:
    945             if type(a) not in types: a = ctx.convert(a)
    946             if type(b) not in types: b = ctx.convert(b)

/usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in <genexpr>(.0)
    252         case the quadrature rule is able to reuse them.
    253         """
--> 254         return self.ctx.fdot((w,w) in nodes)
    255 
    256 

<ipython-input-440-7ebeff3611f6> in <lambda>(U)
      1 def f(n):
----> 2     integrand = lambda U: dxdu_u(U) * mp.besselj(n,nth)
      4     return mp.quadosc(integrand,zeros=bjz)
      5 

这时它开始使用scipy插值代码

/usr/local/lib/python3.6/dist-packages/scipy/interpolate/fitpack2.py in __call__(self,x,nu,ext)
    310             except KeyError:
    311                 raise ValueError("Unknown extrapolation mode %s." % ext)
--> 312         return fitpack.splev(x,self._eval_args,der=nu,ext=ext)
    313 
    314     def get_knots(self):

/usr/local/lib/python3.6/dist-packages/scipy/interpolate/fitpack.py in splev(x,tck,der,ext)
    366         return tck(x,extrapolate=extrapolate)
    367     else:
--> 368         return _impl.splev(x,ext)
    369 
    370 

/usr/local/lib/python3.6/dist-packages/scipy/interpolate/_fitpack_impl.py in splev(x,ext)
    596         shape = x.shape
    597         x = atleast_1d(x).ravel()
--> 598         y,ier = _fitpack._spl_(x,t,c,k,ext)
    599 
    600         if ier == 10:

TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'

_fitpack._spl_可能是编译后的代码（为了提高速度）。它不能直接使用mpmath对象；它必须通过与C兼容的double传递其值。

为说明问题，请制作一个mpmath对象的numpy数组：

In [444]: one,two = mp.mpmathify(1),mp.mpmathify(2)                            
In [445]: arr = np.array([one,two])                                             
In [446]: arr                                                                   
Out[446]: array([mpf('1.0'),mpf('2.0')],dtype=object)

In [447]: arr.astype(float)    # default 'unsafe' casting                                                     
Out[447]: array([1.,2.])
In [448]: arr.astype(float,casting='safe')                                     
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-448-4860036bcca8> in <module>
----> 1 arr.astype(float,casting='safe')

TypeError: Cannot cast array from dtype('O') to dtype('float64') according to the rule 'safe'

使用integrand = lambda U: dxdu_u(float(U)) * mp.besselj(n,U)，

In [453]: f(0)      # a minute or so later                                                                
Out[453]: mpf('0.61060303588231069')