Implement double-exponential radial quadrature about integratorxx HOT 4 OPEN

Great, if you have the closed-form for these quadratures, it should be an easy add.

This also brings up a point I've been meaning to address - we should have stock generators for SG-X spherical quadratures. I'll open an issue to track.

from integratorxx.

This also brings up a point I've been meaning to address - we should have stock generators for SG-X spherical quadratures. I'll open an issue to track.

Yeah, although I think those grids also go hand-in-hand with specific radial quadratures and atomic partitionings. The latter can probably be changed, the former shouldn't. But the whole idea about a standard quadrature library is that these quadratures can be made available to any code.

from integratorxx.

@susilehtola I wanted to hash some of this out before spending time on implementing it. Let's take DE1 as an example

$$ r(x) = \exp( \alpha \sinh(x) ) $$

$$ r'(x) = \alpha \cosh(x)\ r(x) $$

The integral is then approximated (I'm neglecting the $r^2$ spherical Jacobian due to the code conventions)

$$ \int_0^\infty f(r)\mathrm{d}r = \sum_{i=-\infty}^{+\infty} w_i r'(x_i) f(r(x_i)) $$

where $x_i$ is taken to be the uniform trapezoid with grid spacing $h$ (i.e. $w_i = h$). This expression is an infinite sum that must be truncated in practice.

In the paper, they do this by selecting a $R_\mathrm{min}$ and $R_\mathrm{max}$ and truncating the sum to include $x_i\in [r^{-1}(R_\mathrm{min}),\ r^{-1}(R_\mathrm{max})]$. For each of the DE grids, the radial transformation is analytically invertible for $r>0$, e.g.

$$ r^{-1}(R) = \mathrm{asinh}\left(\frac{\ln R}{\alpha}\right) $$

where asinh has an STL implementation.

My understanding of how this should work is the following:

1. Select a sensible $\alpha$, the paper has some defaults, we can make it configurable defaulting to these values
2. Determine the $x$-range via selecting a sensible $R_\mathrm{min}$ and $R_\mathrm{max}$ s.t. $x_\mathrm{min/max} = r^{-1}(R_\mathrm{min/max})$
3. Generate a UniformTrapezoid base quadrature on $x_i \in [x_\mathrm{min}, x_\mathrm{max}]$ with a specified number of points
4. Apply the DE-X transformation to $x_i$ ala the modular design introduced in #27

Make sense?

from integratorxx.

For some undocumented $\alpha$ (but likely $\alpha = \frac{\pi}{2}$ per the paper: edit: confirmed numerically), a select set of nodes and weights are provided in boost::math

from integratorxx.