nintegrate function



On 5/23/07, Stavros Macrakis <macrakis at alum.mit.edu> wrote:
> 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.

I would prefer to have a separate function for numerical integration.

> 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.

This depends on what you want nintegrate to do. I agree that it would
be nicer for nintegrate(0,x,a,b) to return 0. The second one should
return unchanged. If q is a function then nintegrate(q,x,a,b) is
equivalent to nintegrate(q(x),x,a,b). This is very useful for
functions which you can compute for numerical values but q(x) produces
an error.

> 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.

I would like to see a helpful message instead of the error code. Something like

  Numerical integration failed with error code 1 (too many
subintervals were done)
  Try this or that to fix the problem.

>       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?

These are really limitations of the quadpack functions. The first one
happens because the exact value is 0 (it happens for all cases where
the exact value is 0). The second one can be done if you increase the
number of subintervals. Error code 5 says that convergence is too
slow. I don't know how to get maxima to compute it.

> 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;

I can try to work on these examples.

-- 
Andrej