mpmath与scipy的关系：E型

Question

Scipy有很多特殊的函数，特别是贝塞尔函数jn（总是用大写字母J琰n（x）表示）和球形贝塞尔函数spherical_jn（用小写字母J琰n（x）表示）。另一方面，mpmath有quadosc，一种积分快速振荡函数的特殊方法，如jn和spherical_jn。我得到的问题是，似乎mpmath的quadosc不支持，例如scipy的jn作为输入来进行积分。我的意思是，如果我使用从numpy导入的quad，那么就不会得到TypeError的错误，但是quad在x非常大的时候，不太适合计算Jĩn（x）或Jĩn（x）的积分。你知道吗

（***）在通过“振荡求积（quadosc）”找到的站点SymPy中，这个例子来自那里。你知道吗

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.

但如果我用的是从numpy进口的Jïn（x）：

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)

然后我得到了以下错误（编辑：添加了回溯）

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.

那么，如何使用来自mpmath的quadosc中scipy的函数jn？如何修复此错误？谢谢你的帮助。你知道吗

Answer 1

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

Answer 2

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 normal 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 double编译的，但不是更高精度的float。你知道吗

mpmath是否谈到将其与numpy一起使用？我见过mpmath和sympy一起使用，但没有和numpy一起使用。你知道吗