On 5/23/07, Andrej Vodopivec <andrej.vodopivec at gmail.com> wrote:
>
> I would like to see a nice wrapper function for quadpack integration
> methods. I have written a nintegrate function:
Andrej, thanks for your contribution of nintegrate, a simplifying function
which uses numerical methods. I agree that a good wrapper function for the
various numerical routines we have is a worthy goal. I understand this is a
work-in-progress, and I hope you'll find my comments below constructive.
A big issue in the design is: should nintegrate in fact be a separate
function from integrate? Perhaps integrate should try numerical techniques
if symbolic ones fail, under control of a flag (e.g. numer)? This would
generalize nicely to other parts of the system, e.g. numeric limits, numeric
sums, etc.
The current nintegrate in fact doesn't simplify at all if the limits aren't
numeric, so you get
nintegrate(0,x,a,b) => no change
nintegrate(q,x,a,b) => no change
cases which the normal nounform integrate simplifies.
I assume you plan to improve the error returns to give a
human-understandable error rather than a number (and to suppress the
Fortran-style cruft -- it would be nice to suppress the compilation cruft as
well), and in the cases where it is a limitation in the underlying routines
rather than a user error, give back the nintegrate expression, e.g.
nintegrate(signum(x),x,-1,1) => error code 2
(why is this failing? discontinuity??)
nintegrate(sin(1/x),x,0.00001,1) => error code 1
(OK, this is a hard one, but x*sin(1/x) and abs(sin(1/x)) also
fail)
nintegrate(x^(-10001/10000),x,1,inf) => error code 5
(symbolic integration can do this one easily)
I don't know much about numeric integration routines; should nintegrate be
using other, more robust ones in some of the above cases?
There are also some bugs in the current version that need to be fixed:
nintegrate(f(x),x,1,2) => internal error
should just return with no change
nintegrate(y,x,1,2) => undefined function y
should return either unchanged, or y
and should never try to interpret y as a function,
which is inconsistent with integrate
nintegrate(y*x,x,1,2) => internal error, should return as is
nintegrate(1,x,0,2^%e) => returns as is; should bind $numer:true
when $float'ing
(arguably a bug in float)
I hope this is helpful.
-s