Comparison to quaddeo?

[komatsu5147]

Very nice! I wonder how this method compares with quaddeo of https://github.com/machakann/DoubleExponentialFormulas.jl Any chance to make the comparison?

[xzackli]

Hey @komatsu5147,

I'm also curious about this, and I plan to make some work-precision plots for this package next week. I'll ping you when I have them. It should be trivial to derive the relation for sine and cosine transforms needed for Levin methods, to make a "strategy" comparison with quaddeo. One important difference is that although DE can be extended to handle Bessel integrals, I don't believe DoubleExponentialFormulas.jl currently supports them.

The Levin-type method of this package is intended for finite 0 < a < b with no singularities, so DE may be a better choice when such conditions do not hold. However, this Levin-type method offers a very flexible integration mesh, which one can manually enforce by passing timesteps to the ODE solver as a keyword argument. I wrote this package so that Bolt can have a high-accuracy method for computing the line-of-sight integral of the photon source function. This conformal time integral contains a rapidly-oscillating spherical Bessel function multiplied by the photon perturbation. We need high-resolution quadrature at last scattering (huge contribution in a narrow shell), but we can afford 10x lower resolutions after recombination finishes.

[komatsu5147]

Thanks! I think that quaddeo can handle the Bessel integrals because all Bessel functions asymptote to a sin(x+a) function with some phase shift a when the argument x is large, which is all you need to apply quaddeo. Eager to know the comparison, especially the computational cost. quaddeo is powerful but it tends to require many evaluations. Perhaps the Levin method requires less. I did not know about Bolt.jl - will check it out!

[xzackli]

Hey @komatsu5147, I was recently reminded of the series of papers Piessens & Branders 1983, a simple Chebysev-based quadrature. I've seen this strategy used in cosmology recently, in Campagne+17 (AngPow). Cosmologists tend to do the same integral many times since we like to sample, so this quadrature method can be very competitive if one precomputes and stores the weights. If one uses the same grid every time, then one only evaluates the function on the Chebyshev nodes of each subinterval. My "bessel integrator race" will also involve a quick implementation of Piessens & Branders.

I am trying to persuade some people I know to write a paper along these lines. I think the Chebyshev basis is under-utilized, as it has some very nice numerical properties.

[komatsu5147]

Interesting! I've used the Chebyshev basis for interpolation of a function, which then allows for analytical integration and differentiation, but have not used it in combination with the oscillatory integrals. Thanks for the references.