2

我需要对 g(u)jn(u) 类型进行积分，其中 g(u) 是一个没有零的平滑函数，而贝塞尔函数中的 jn(u) 具有无穷大的零，但出现以下错误：

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

首先，我需要将变量 x 更改为变量 u 并在新变量 u 中进行积分，但是函数 u(x) 如何在解析上不可逆，因此我需要使用插值来进行数值求逆。

import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline

x = np.linspace(0.1, 100, 1000)
u = lambda x: x*np.exp(x)
dxdu_x = lambda x: 1/((1+x) * np.exp(x))               ## dxdu as function of x: not invertible
dxdu_u = InterpolatedUnivariateSpline(u(x), dxdu_x(x)) ## dxdu as function of u: change of variable

在此之后，积分为：

from mpmath import mp

def f(n):
    integrand = lambda U: dxdu_u(U) * mp.besselj(n,U)
    bjz = lambda nth: mp.besseljzero(n, nth)
    return mp.quadosc(integrand, [0,mp.inf], zeros=bjz)

我使用quadoscfrommpmath而不是quadfromscipy因为quadosc更适合于对快速振荡函数进行积分，例如 Bessel 函数。但是，另一方面，这迫使我使用两个不同的包，通过插值计算，scipy计算贝塞尔函数和产品的积分，所以我怀疑两个不同包的混合会产生一些问题/冲突。所以当我做：dxdu_umpmathmp.besselj(n,U)dxdu_u(U) * mp.bessel(n,U)

print(f(0))

我得到了错误：

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-38-ac2976a6b736> in <module>
     12     return mp.quadosc(integrand, [0,mp.inf], zeros=bjz)
     13 
---> 14 f(0)

<ipython-input-38-ac2976a6b736> in f(n)
     10     integrand = lambda U: dxdu_u(U) * mp.besselj(n,U)
     11     bjz = lambda nth: mp.besseljzero(n, nth)
---> 12     return mp.quadosc(integrand, [0,mp.inf], zeros=bjz)
     13 
     14 f(0)

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

有谁知道我该如何解决这个问题？谢谢

4

2 回答 2

2

完整的回溯（您狙击的部分）表明错误出__call__在 univariatespline 对象的方法中。所以确实问题在于 mpmath 集成例程输入mpf小数，而 scipy 无法处理它们。

一个最简单的解决方法是手动将参数的违规部分integrand转换为浮点数：

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

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

一种可能的替代方法可能是避免 mpmath 并使用weights参数 to scipy.integrate.quad，请参阅文档（向下滚动到weights="sin"部分）

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

于 2020-01-08T20:15:33.747 回答
0

完整的追溯：

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, f, interval, omega=None, period=None, zeros=None):

/usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in quad(ctx, f, *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, f, points, prec, epsilon, 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, degree, prec, results, verbose))
    233                 if degree > 1:
    234                     err = self.estimate_error(results, prec, epsilon)

/usr/local/lib/python3.6/dist-packages/mpmath/calculus/quadrature.py in sum_next(self, f, nodes, degree, prec, 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, f(x)) for (x,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,U)
      3     bjz = lambda nth: mp.besseljzero(n, nth)
      4     return mp.quadosc(integrand, [0,mp.inf], 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, der, extrapolate=extrapolate)
    367     else:
--> 368         return _impl.splev(x, tck, der, ext)
    369 
    370 

/usr/local/lib/python3.6/dist-packages/scipy/interpolate/_fitpack_impl.py in splev(x, tck, der, ext)
    596         shape = x.shape
    597         x = atleast_1d(x).ravel()
--> 598         y, ier = _fitpack._spl_(x, der, 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 兼容的双精度值传递。

为了说明这个问题，制作一个 numpympmath对象数组：

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')
于 2020-01-08T20:57:49.810 回答