Связь между mpmath и scipy: ошибка типа

В Scipy есть много специальных функций, в частности функции Бесселя jn (всегда обозначаются прописной буквой J_n(x)) и сферические функции Бесселя spherical_jn (обозначаются строчной буквой j_n(x)). С другой стороны, в mpmath есть quadosc, специальный метод для интеграции быстро колеблющихся функций, таких как jn и spherical_jn. Проблема, которую я обнаружил, заключается в том, что кажется, что quadosc из mpmath не поддерживается, например, jn из scipy в качестве входных данных для создания этого интеграла. Я имею в виду, что если я использую quad, импортированный из numpy, поэтому не возникает никакой ошибки TypeError, но quad не очень подходит для оценки интегралов J_n (x) или j_n (x), когда x очень большой.

(***) На сайте SymPy найти по "Колебательной квадратуре (quadosc)", этот пример взят оттуда.

from mpmath import findroot, quadosc, inf, j0

j0zero = lambda n: findroot(j0, pi*(n-0.25)) # ***
I = quadosc(j0, [0, inf], zeros=j0zero)
print(I)
I = 1.0 # OK, this is the correct answer.

Но если я использую J_n(x), импортированный из numpy:

from scipy.special import jn

f = lambda x: jn(0,x)
j0zero = lambda n: findroot(f, pi*(n-0.25))
II = quadosc(f, [0, inf], zeros=j0zero)
print(II)

затем я получил следующую ошибку (отредактировано: добавлен Traceback)

TypeError Traceback (most recent call last)
~/anaconda3/lib/python3.7/site-packages/mpmath/calculus/optimization.py in findroot(ctx, f, x0, solver, tol, verbose, verify, **kwargs)
    927  try:
--> 928   fx = f(*x0)
    929    multidimensional = isinstance(fx, (list, tuple, ctx.matrix))
<ipython-input-449-aeebd9a1e908> in <lambda>(x)
      2 
----> 3 f = lambda x: jn(0,x)
      4 j0zero = lambda n: findroot(f, pi*(n-0.25))

TypeError: ufunc 'jv' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

During handling of the above exception, another exception occurred:
TypeError  Traceback (most recent call last)
<ipython-input-449-aeebd9a1e908> in <module>
      3 f = lambda x: jn(0,x)
      4 j0zero = lambda n: findroot(f, pi*(n-0.25))
----> 5 II = quadosc(f, [0, inf], zeros=j0zero)
      6 print(II)

~/anaconda3/lib/python3.7/site-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)]

~/anaconda3/lib/python3.7/site-packages/mpmath/calculus/optimization.py in findroot(ctx, f, x0, solver, tol, verbose, verify, **kwargs)
    929     multidimensional = isinstance(fx, (list, tuple, ctx.matrix))
    930     except TypeError:
--> 931     fx = f(x0[0])
    932     multidimensional = False
    933     if 'multidimensional' in kwargs:

С другой стороны, если я использую quad, то я получаю

from scipy.integrate import quad

f = lambda x: jn(0,x)
III = quad(f,0,inf)[0]
print(III)
III = -21.154674722694516 # What is an incorrect answer.

Итак, как я могу использовать функцию jn из scipy внутри quadosc из mpmath? Как я могу исправить эту ошибку? Спасибо за помощь.


python numpy scipy mpmath
person Diogo    schedule 22.12.2019    source источник
comment
Поскольку вы не используете jv напрямую, покажите трассировку, чтобы мы могли увидеть, где возникает ошибка. Я подозреваю, что mpmath передает свой собственный класс чисел функции scioy.numpy, которую он не может обработать.   -  person hpaulj    schedule 22.12.2019
comment
Спасибо, я добавил трассировку. Да, мой вопрос в том, что существует какой-то способ заставить mpmath понимать scipy-функции.   -  person Diogo    schedule 22.12.2019
comment
Важно сказать, что я тестировал jn и jv функцию scipy для случая n=v=0, то есть jn(0,x) и jv(0,x) и обе они являются одной и той же функцией.   -  person Diogo    schedule 22.12.2019
comment
да, глядя на документы jn, я вижу, что это другое имя для jv.   -  person hpaulj    schedule 22.12.2019
comment
Сначала я хотел понять, почему ошибка говорит jv, когда вы звоните jn. Это прояснилось. Проблема в том, что quadosc передает mpmath числа функции jn, которая может использовать только C удвоений (или эквивалентных чисел с плавающей запятой).   -  person hpaulj    schedule 22.12.2019
comment
Да, но проблема не в этом, даже если я поменяю jn на jv, проблема не исчезнет, ​​и я получу ту же ошибку. Проблема в том, что quadosc не понимает jn (или jv) как интегрант. Таким образом, если существует какой-то способ заставить mpmath и scipy взаимодействовать между собой, то quadosc можно применять для большого класса специальных функций, особенно осциллирующих функций.   -  person Diogo    schedule 22.12.2019

Ответы (2)


arrow_upward
1
arrow_downward 

In [31]: one,two = mpmath.mpmathify(1), mpmath.mpmathify(2)                     

In [32]: one,two                                                                
Out[32]: (mpf('1.0'), mpf('2.0'))

In [33]: one+two                                                                
Out[33]: mpf('3.0')

In [34]: jn(1,2)                                                                
Out[34]: 0.5767248077568736

In [35]: jn(one,two)                                                            
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-35-ec48c25f686b> in <module>
----> 1 jn(one,two)

TypeError: ufunc 'jv' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

mpmath числа нельзя использовать в jn/jv.

Точно так же число mpmath может использоваться в функции mpmath, но не эквивалентно numpy:

In [41]: mpmath.sin(one)                                                        
Out[41]: mpf('0.8414709848078965')

In [42]: np.sin(1)                                                              
Out[42]: 0.8414709848078965

In [43]: np.sin(one)                                                            
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
AttributeError: 'mpf' object has no attribute 'sin'

The above exception was the direct cause of the following exception:

TypeError                                 Traceback (most recent call last)
<ipython-input-43-38fc918b311a> in <module>
----> 1 np.sin(one)

TypeError: loop of ufunc does not support argument 0 of type mpf which has no callable sin method

Вы можете преобразовать обычный номер Python mpmath с помощью:

In [44]: float(one)                                                             
Out[44]: 1.0

In [45]: jn(float(one),float(two))                                              
Out[45]: 0.5767248077568736

np.float64(one) тоже работает, но jn не нравится np.float128(one). Очевидно, что jn был скомпилирован для C doubles, но не для float с более высокой точностью.

Говорит ли mpmath об использовании его с numpy? Я видел, как mpmath использовался с sympy, но не с numpy.

person hpaulj    schedule 22.12.2019

arrow_upward
0
arrow_downward

Невозможно использовать scipy.special.jv без преобразования его аргументов в числа с плавающей запятой/целые числа.

person ev-br    schedule 23.12.2019