Open haasal opened 1 day ago
Thank you for this excellent and well-thought-out proposal! You raise several compelling points, especially regarding the limitations of current integration methods. You're absolutely right about the challenges with Dirac delta-like distributions - most other integration techniques indeed fail in these cases, and the current implementation does struggle with both accuracy and performance over large intervals.
The benefits of QAGS sound quite promising based on your description. While I haven't had the chance to look into it in detail yet, the ability to handle singularities effectively is a significant advantage. Additionally, your suggestion to implement QAGI for infinite interval integration would fill an important gap in our current functionality, as we currently don't have a robust solution for this case.
However, I should be transparent about the timeline: due to current workload, I likely won't be able to start implementing this until December. That said, if you're interested in implementing this yourself, we would definitely welcome a PR! Your understanding of both the mathematical background and the practical benefits of QAGS/QAGI would be valuable in this implementation.
Thanks again for bringing this up and offering to help with the implementation. It's exactly these kinds of well-reasoned proposals that help improve the library's functionality and performance.
Very nicely written response ;)
I would use the gsl implementations as a reference. The python c-implementation is horrible to read. But pythons QAGS implemtentation is different from gsl's implementation. For example in gsl QAGS I can not integrate from -inf to inf (only using QAGI but this resulted for me in oscillations probably due to numerical errors). In scipy you can enter -inf to inf as integration boundries into quad(...)
.
This whole thing hast cost me countless hours over the last three days and I still haven't found a solution to my particular problem.
I'm trying to fit superconducting gaps which are notoriously difficult so maybe this is the issue.
My mathematical background in numerical integration isn't particularly versed but I will probably try to read up on these algorithms.
Reference implementation.
Note: In scipy they have default handlers for non-convergence or round-off errors and in these cases the integrator simply returns nan
(I think). I propose simply returning a Result instead.
I propose following and extending scipys quad function which is based on QUADPACKS adaptive QAGS algorithm. I also propose setting the default integrate function from the prelude to this one.
The reason behind this is as follows:
In short: there is a reason the scipy project decided on QAGS as their quadrature algo of choice.
I would take the implementation of the gsl project.
This probably wouldn't be too hard to implement as the original algorithm is based on GK which we already have.
I would of course help if this gets accepted!