python - Relation between mpmath and scipy: Type Error

Question

Scipy have a lot of special functions, in particular Bessel functions jn (always denoted by uppercase letter J_n(x)) and spherical Bessel functions spherical_jn (denoted by lowercase letter j_n(x)). On the other hand mpmath have quadosc, a special method for integrate rapidly oscillating functions, like jn and spherical_jn. The problem I had obtained is that seems that quadosc from mpmath not support, e.g, jn from scipy as an input to make this integral. I means, if I use quad imported from numpy, so nothing error of TypeError is obtained, but quad is not very appropriate to evaluate integrals of J_n(x) or j_n(x) when x is very large.

(***) In the site SymPy find by "Oscillatory quadrature (quadosc)", this example come from there.

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.

But if I use J_n(x) imported from 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)

then I got the following error (Edited: added the 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:

On the other hand, if I use quad then I got

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.

So how can I use a function jn that come from scipy inside a quadosc of mpmath? How can I fix this error? Thanks for all help.

Answer 1

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

您可以使用以下方式转换 mpmath 普通 python 数字：

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 双精度数编译，但不是更高精度的浮点数。

是否mpmath谈论使用它numpy？我见过mpmath使用 with sympy，但不是 with numpy。

Answer 2

scipy.special.jv如果不将其参数转换为浮点数/整数，就无法使用。