TypeError: Невозможно привести данные массива из dtype ('O') к dtype ('float64') в соответствии с правилом safe.

Question

Мне нужно сделать интеграл типа 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)

Я использую quadosc от mpmath, а не quad от scipy, потому что 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'

Кто-нибудь знает, как я могу решить эту проблему? Спасибо

Answer 1

Полный возврат (часть, которую вы отсканировали) показывает, что ошибка в методе __call__ объекта univariatespline. Таким образом, проблема в том, что подпрограмма интеграции mpmath подает свои mpf десятичные дроби, и scipy никак не может с ними справиться.

Самое простое исправление заключается в том, чтобы вручную приводить вызывающую ошибку часть аргумента integrand до числа с плавающей запятой:

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

Как правило, это подвержено числовым ошибкам (mpmath специально использует высокоточные переменные!), Поэтому действуйте с осторожностью. В этом конкретном случае c все может быть в порядке, потому что интерполяция фактически выполняется с двойной точностью. Тем не менее, лучше всего проверить результаты.

Возможной альтернативой может быть предотвращение mpmath и использование аргумента weights для scipy.integrate.quad, см. документы (прокрутите вниз до части weights="sin" )

Другая альтернатива - полностью придерживаться mpmath и реализовать интерполяцию в чистом python (таким образом, mpf объекты, вероятно, хороши, поскольку они должны поддерживать обычную арифметику). Вероятно, достаточно простой линейной интерполяции. Если это не так, то не так уж сложно написать свой собственный кубический c сплайн-интерполятор.

Answer 2

Полный возврат:

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 совместимых двойников.

Чтобы проиллюстрировать проблему, создайте массив numpy из mpmath объектов:

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')