类型错误:无法根据规则“安全"将数组数据从 dtype('O') 转换为 dtype('float64') [英] TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'

问题描述

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

I need to make an integral of the type g(u)jn(u) where g(u) is a smooth function without zeros and jn(u) in the Bessel function with infinity zeros, but I got the following error:

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

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

First I need to change of variable x to variable u and make an integration in the new variable u but how the function u(x) is not analytically invertible so I need to use interpolation to make this inversion numerically.

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)

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

I use quadosc from mpmath and not quad from scipy because quadosc is more appropriate to make integral of rapidly oscillating functions, like Bessel functions. But, by other hand, this force me to use two different packges, scipy to calculate dxdu_u by interpolation, and mpmath to calculate the Bessel functions mp.besselj(n,U) and the integral of the product dxdu_u(U) * mp.bessel(n,U) so I suspect that this mix of two different packages can make some issue/ conflict. So when I make:

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'

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

Does anyone know how I can solve this problem? Thanks

推荐答案

完整的回溯(你狙击的部分)显示错误在 univariatespline 对象的 __call__ 方法中.所以确实问题在于 mpmath 集成例程以它的 mpf 小数形式输入，而 scipy 无法处理它们.

The full traceback (the part you sniped) shows that the error is in the __call__ method of the univariatespline object. So indeed the problem is that the mpmath integration routine feeds in its mpf decimals, and scipy has no way of dealing with them.

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

A simplest fix is then to manually cast the offending part of the argument of the integrand to a float:

被积函数 = lambda U: dxdu_u(float(U)) * mp.besselj(n,U)

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

In general this is prone to numerical errors (mpmath uses its high-precision variables on purpose!) so proceed with caution. In this specific case it might be OK, because the interpolation is actually done in double precision. Still, best check the results.

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

A possible alternative might be to avoid mpmath and use the weights argument to scipy.integrate.quad, see the docs (scroll down to weights="sin" part)

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

Another alternative is to stick with mpmath all the way and implement the interpolation yourselves in pure python (this way, mpf objects are probably fine since they should support usual arithmetics). It's likely a simple linear interpolation is enough. If it's not, it's not too big of a deal to code up your own cubic spline interpolator.

这篇关于类型错误:无法根据规则“安全"将数组数据从 dtype('O') 转换为 dtype('float64')的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！