1

I have a density function (from quantum mechanics calculations) to be multiplied with the spherical Bessel function with a momentum grid (momentum q 1d array, real space distance r 1d array, so need to calculate jn(q*r) 2d array). The product will be integrated across the real space to get results as function of momentum (results in 1d array same shape as q).

The Bessel function has oscillation, while the density function fast decay over a threshold distance. I used the adaptive quadrature in quadpy which is fine when oscillation is slow but it fails with high oscillation (for high momentum values or high orders in Bessel functions). The mpmath quadosc could be a nice option, but currently I have the problem that "object arrays are not supported", which seems to be the same as in Relation between mpmath and scipy: Type Error, what would be the best way to solve it since the density function is calculated outside of the mpmath.

from mpmath import besselj, sqrt, pi, besseljzero, inf,quadosc
from scipy.interpolate import interp1d

n = 1
q = np.geomspace(1e-7, 500, 1000)
# lets create a fake gaussian density 
x = np.geomspace(1e-7, 10, 1000)
y = np.exp(-(x-5)**2)
density = interp1d(x,y,kind='cubic',fill_value=0,bounds_error=False)


# if we just want to integrate the spherical bessel function
def spherical_jn(x,n=n):
    return besselj(n + 1 / 2, x) * sqrt(pi / 2 / x)
# this is fine
vals = quadosc(
    spherical_jn, [0, inf], zeros=lambda m: besseljzero(n + 1 / 2, m)
)

# now we want to integrate the spherical bessel function times the density
def spherical_jn_density(x,n=lprimeprime):
    grid = q[..., None] *x
    return besselj(n + 1 / 2,  grid) * sqrt(pi / 2 / grid)*density(x)

# this will fail
vals_density = quadosc(
    spherical_jn_density, [0, inf], zeros=lambda m: besseljzero(n + 1 / 2, m)
)

Expect: an accurate integral of highly oscillating spherical Bessel function with arbitrary decaying function (decay to zero at large distance).

Zezhong Zhang
  • 13
  • 3
  • Most (all) `scipy` is based on `numpy`, and the compiled parts use the C numeric types - float, double. An array of `mpmath` numbers will be object dtype, the elements being those `mpmath` objects, not more limited precision `doubles`. So you have to `mpmath` code through out, and not depend on anything from `scipy`. – hpaulj Feb 06 '23 at 17:04
  • The density itself is not calculated with scipy but external program that returns arrays. Here I interface it with numpy and make scipy interp1d so it can act as a function. Is it possible to convert my density matrix (a 1d numpy array as function of distance) and grid (as function of momentum and distance, 2d numpy array) to a mpmath function? If not, I think then I have to give up the idea of using mpmath... – Zezhong Zhang Feb 06 '23 at 18:58
  • I probably should have asked for a full error traceback, rather than try to explain the imagined incompatibility. But I will guess that `densiry(x)` is causing the problem when given a `mpmath` object. READ YOUR ERROR MESSAGE CAREFULLY – hpaulj Feb 06 '23 at 19:07

1 Answers1

2

Your density is interp callable, which works like:

In [33]: density(.5)
Out[33]: array(1.60522789e-09)

It does not work when given a mpmath object:

In [34]: density(mpmath.mpf(.5))
ValueError: object arrays are not supported

It's ok if x is first converted to ordinary float:

In [37]: density(float(mpmath.mpf(.5)))
Out[37]: array(1.60522789e-09)

Tweaking your function:

def spherical_jn_density(x,n=1):
    print(repr(x))
    grid = q[..., None] *x
    return besselj(n + 1 / 2,  grid) * sqrt(pi / 2 /grid) * density(x)

and trying to run the quadosc (with a smaller q)

In [57]: vals_density = quadosc(
    ...:     spherical_jn_density, [0, inf], zeros=lambda m: besseljzero(n + 1 / 2, m))
mpf('0.506414729137261838698106')
TypeError: cannot create mpf from array([[mpf('5.06414729137261815781894e-8')],
       [mpf('0.000000473559111442409924364745')],
       [mpf('0.00000442835129247081824275722')],
       [mpf('0.0000414104484439061558283487')],
       [mpf('0.000387237851532012775822723')],
       [mpf('0.00362114295531604773233197')],
       [mpf('0.0338620727569835882851491')],
       [mpf('0.316651395857188250996884')],
       [mpf('2.96107409661232278850947')],
       [mpf('27.6896294168963266721213')]], dtype=object)

In other words,

besselj(n + 1 / 2,  grid)

is having problems, even before trying to evaluate density(x). mpmath functions don't work with numpy arrays; and many numpy/scipy functions don't work with mpmath objects.

hpaulj
  • 221,503
  • 14
  • 230
  • 353