0

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.

hpaulj
  • 221,503
  • 14
  • 230
  • 353
Diogo
  • 64
  • 2
  • 7
  • Since you don't use `jv` directly, show the traceback so we can see where the error occurs. I suspect mpmath is passing its own class of numbers to the scioy.numpy function, which it cannot handle. – hpaulj Dec 22 '19 at 02:38
  • Thanks, I added the traceback. Yes, my question is there exist some way of make mpmath understand scipy functions. – Diogo Dec 22 '19 at 03:49
  • It is important to say that I tested `jn` and `jv` scipy function for the case n=v=0, that is, `jn(0,x)` and `jv(0,x)` and both are the same function. – Diogo Dec 22 '19 at 04:04
  • yes, looking at the `jn` docs, I see that that's another name for `jv`. – hpaulj Dec 22 '19 at 04:27
  • Initially I wanted to understand why the error said `jv` when you called `jn`. That's cleared up. The problem is that `quadosc` passes `mpmath` numbers to the `jn` function, which can only use `C` doubles (or equivalent floats). – hpaulj Dec 22 '19 at 04:43
  • Yes, but this is not the problem, even if I change of `jn` to `jv` the problem persist and I got the identical error. The problem is that quadosc not understand `jn` (or `jv`) as an integrand. So if there exist some way of make mpmath and scipy communicate between, then quadosc can be applied for a large class of special functions, specially oscillating functions. – Diogo Dec 22 '19 at 04:46

2 Answers2

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 numbers cannot be use in jn/jv.

In the same way, a mpmath number can be used in a mpmath function, but not the equivalent 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

You can convert the mpmath normal python number with:

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

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

np.float64(one) also works, but jn does not like np.float128(one). Evidently jn has been compiled for C doubles, but not higher precision floats.

Does mpmath does talk about using it with numpy? I've seen mpmath used with sympy, but not with numpy.

hpaulj
  • 221,503
  • 14
  • 230
  • 353
0

There is no way to use scipy.special.jv without converting its arguments to floats/integers.

ev-br
  • 24,968
  • 9
  • 65
  • 78