f0uriest / quadax

Numerical quadrature with JAX
MIT License
39 stars 1 forks source link

Question about absolute error in NestedRule #13

Open e-eight opened 1 month ago

e-eight commented 1 month ago

Nice package! I am trying to implement some custom quadrature rules using quadax. I was trying to understand how you are computing the absolute error in NestedRule:

abserr = jnp.where(
    (integral_mmn != 0.0) & (abserr != 0.0),
    integral_mmn * jnp.minimum(1.0, (200.0 * abserr / integral_mmn) ** 1.5),
    abserr,
)
abserr = jnp.where(
    (integral_abs > uflow / (50.0 * eps)),
    jnp.maximum((eps * 50.0) * integral_abs, abserr),
    abserr,
)

Where are you getting the numbers 200.0, 1.5, and 50.0 from?

f0uriest commented 1 month ago

These were mostly taken from the fortran library QUADPACK (for which the documentation is rather sparse for things like this). My understanding is that they are heuristic values meant to compensate for floating point errors and possible under/overflow.

I think for most reasonably "well behaved" integrands this sort of stuff isn't necessary, so you can probably omit it.