double exponential integrator
I'm soliciting feedback on an implementation of integration using the double exponential transform. I've tentatively placed it in scipy.integrate as de_integrate. It's available at this url and branch: http://broad.mit.edu/~thouis/scipy.git DEintegrator (assuming I've set up my git repository correctly, which is quite possibly not the case.) See http://www.johndcook.com/double_exponential_integration.html for an explanation of the double exponential transform. It's main use is integrating functions that have a singularity at one or both ends of the integration limits, though it works reasonably well for arbitrary functions, though my implementation could be made more robust in many ways. I've also included code for generating the constants necessary for integrations with infinite limits (0 to infinity as well as -inf to +inf). Thouis Jones
Hi, Thu, 19 Mar 2009 13:55:58 -0400, Thouis (Ray) Jones wrote:
I'm soliciting feedback on an implementation of integration using the double exponential transform. I've tentatively placed it in scipy.integrate as de_integrate.
+1 for including this in Scipy 0.8.0. Some quick comments (I'll try to find time for better comments later): - Vectorization: you are using a list comprehension in doubleexp.py:contribution_at_level and elsewhere These statements could be vectorized -- I believe you can also require that the integrand function `f` can evaluate many points at the same time and return an array. Could be an useful speedup. - Is generating the _abscissas_raw and _weights_raw costly? I see that you use `mpmath` to prepare these. Does the generation fail in double precision? If not, it might be better to generate them when the integration function is first used. (Also, it might be nice to put also the generator functions in the same doubleexp.py; these are not long files.) - Function name `de_integrate`: Perhaps it should be `quad_de` or something similar, since it's a quadrature, and the "basic" quadrature in scipy.integrate is called `quad`. - I'm not sure about if a `full_output` switch is good API. Several Scipy functions do currently use something like that, but perhaps it would be better not to introduce more...
It's available at this url and branch: http://broad.mit.edu/~thouis/scipy.git DEintegrator
(assuming I've set up my git repository correctly, which is quite possibly not the case.)
The repository works OK. -- Pauli Virtanen
Renamed the function to quad_de. Changed its signature to be closer to quadrature(). Improved reporting when the func() returns nonfinite values. Added vectorization (via the vectorize1 functions in quadrature.py. Filed a trac ticket: #895. The generation of the weights and abscissas is not particularly costly, but I believe they do fail in double precision arithmetic. If it's really desirable to have them created dynamically (to allow more levels of integration, for instance), it would probably require writing a Taylor expansion or similar. I've left them in a separate file for now, due to their reliance on mpmath. In the future, this code could accept arbitrary limits (0, 1, or 2 infinite limits) and use whichever abscissa and weight functions are appropriate, with C code to generate them on demand to arbitrary levels. I would like to gauge interest in this code in general before adding more functionality. The interface would be backwards compatible, with possibly one new keyword argument to handle the two different versions of integrals from 0 to infinity (standard versus exponential falloff). Ray Jones On Fri, Mar 20, 2009 at 05:46, Pauli Virtanen <pav@iki.fi> wrote:
Hi,
Thu, 19 Mar 2009 13:55:58 -0400, Thouis (Ray) Jones wrote:
I'm soliciting feedback on an implementation of integration using the double exponential transform. I've tentatively placed it in scipy.integrate as de_integrate.
+1 for including this in Scipy 0.8.0.
Some quick comments (I'll try to find time for better comments later):
- Vectorization: you are using a list comprehension in doubleexp.py:contribution_at_level and elsewhere
These statements could be vectorized -- I believe you can also require that the integrand function `f` can evaluate many points at the same time and return an array. Could be an useful speedup.
- Is generating the _abscissas_raw and _weights_raw costly?
I see that you use `mpmath` to prepare these. Does the generation fail in double precision?
If not, it might be better to generate them when the integration function is first used. (Also, it might be nice to put also the generator functions in the same doubleexp.py; these are not long files.)
- Function name `de_integrate`: Perhaps it should be `quad_de` or something similar, since it's a quadrature, and the "basic" quadrature in scipy.integrate is called `quad`.
- I'm not sure about if a `full_output` switch is good API. Several Scipy functions do currently use something like that, but perhaps it would be better not to introduce more...
It's available at this url and branch: http://broad.mit.edu/~thouis/scipy.git DEintegrator
(assuming I've set up my git repository correctly, which is quite possibly not the case.)
The repository works OK.
-- Pauli Virtanen
_______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
Thu, 19 Mar 2009 13:55:58 -0400, Thouis (Ray) Jones wrote:
I'm soliciting feedback on an implementation of integration using the double exponential transform. I've tentatively placed it in scipy.integrate as de_integrate.
Also, could you file an enhancement ticket on Scipy's Trac: http://projects.scipy.org/scipy/ This should ensure that your contribution does not disappear in the mailing list traffic... -- Pauli Virtanen
participants (2)
-
Pauli Virtanen
-
Thouis (Ray) Jones