【发布时间】:2020-01-08 19:39:39
【问题描述】:
我需要对 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)
我使用来自mpmath 的quadosc 而不是来自scipy 的quad,因为quadosc 更适合于对快速振荡函数(如贝塞尔函数)进行积分。但是,另一方面,这迫使我使用两个不同的包,scipy 通过插值计算dxdu_u,mpmath 计算贝塞尔函数mp.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'
有谁知道我该如何解决这个问题? 谢谢
【问题讨论】:
-
dtype('O') 表示您有一个对象数组。我认为 mp.quadosc 返回某种形式的多精度类型,这是 numpy 无法识别的
-
是的,我想这就是重点。我不知道是否存在某种方式使 mpmath 和 numyp 相互通信。
-
当
numpy/scipy函数使用编译代码时,数字输入必须与C 兼容doubles。mpmath对象不是双精度对象。numpy不“知道”任何关于mpmath的信息,所以你真的无法让它们交流,除非将mpmath对象转换为常规数字。 -
@hpulj 确实,很好的发现!不过这里的问题更简单,因为 mpf/double 问题发生在插值中
标签: python scipy typeerror mpmath