search

f0uriest / quadax

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

Weighted quad: algebraic-logarithmic weight functions #8

Closed twhentschel closed 4 months ago

twhentschel commented 6 months ago

This is a still a work in progress, but I think it's pretty close so I wanted to share a draft for now and hopefully get some ideas on how to fix the remaining bugs.

Current issues:

  1. In the adaptive_quadrature method, we map the integration interval into [-1, 1]. When we need to check if we should apply the weighted quadrature rule instead of the default rule, the code checks if the singular point is in the interval or not. The problem is that the singular point has not been mapped to [-1, 1], so this check usually fails or doesn't work right. When I remove the mapping, things seem to be working (for the most part, see next issue).
  2. There have been some simple test cases that are not passing because I'm getting a nan for the error. This is because integral_mmn in

https://github.com/f0uriest/quadax/blob/c4544fbc6c8be47459489e2995a2098437944c4e/quadax/fixed_order.py#L123-L125

can (incorrectly) be negative, so that when it's taken to the 1.5 power in

https://github.com/f0uriest/quadax/blob/c4544fbc6c8be47459489e2995a2098437944c4e/quadax/fixed_order.py#L132-L136

it gives nan.

I think integral_mmn becomes negative because taking the absolute value of the integrand ruins the scheme I'm using for weighted quadrature. When we have weight functions present in addition to $f(x)$, self.wh contains the influence of the weight function. If we want the integral of the absolute value of the whole integrand (as we do in integral_mmn), then we need to consider the absolute value of the weight functions too. But this becomes a "different" weight function and so self.wh is incorrect here (and I'm not sure there is a simple way to compute wh in this case).

Other small things

twhentschel commented 5 months ago

Hey sorry it's been a while, I've been in the process of moving and starting a new job. @f0uriest I was hoping you might have some thoughts on this PR! In particular, if you have any pointers on understanding where the quadrature error estimate comes from, it might help me come up with a way to compute the error when using weight functions. For the first issue, scipy doesn't allow infinite bounds for this type of weighted quadrature, so maybe we could follow them. An easy way out is to only apply the mapping to the interval (-1, 1) if either of the bounds are infinite, and throw an error if we want to use certain weight function in this case. Otherwise just use the provided intervals.

f0uriest commented 4 months ago

No worries, I've also been dealing with moving and starting a new job :)

The abstract rule stuff is now merged into main (you'll need to change the base branch before reopening the PR, it seems I can't do it from my end, and the PR auto-closed when the old base branch was merged/deleted).

I think it's probably fine to only allow weight functions for finite intervals like scipy. I think in theory it should be possible to get it to work correctly for infinite intervals but we can figure that out later.

For the most part the error estimate comes from using two different methods to estimate the integral: one low order rule, and another high order rule. For efficiency all the ones I've implemented are so called "nested rules" where the low and high order rules share the same nodes so function values can be re-used. The difference between the low and high order estimates gives an estimate of the error, which is then used in the outer adaptive_quadrature loop to decide which interval to sub-divide.

The other parts integral_abs and integral_mmn are mostly there to account for the possibility of roundoff error (ie, if the error estimate is less than machine precision it can't really be trusted). The particulars of the code are basically taken from quadpack, I haven't found any rigorous justification or explanation but they seem to work well in practice.