bug in integration

[ggrafendorfer]

Hi,

sage 5.0 beta11, on Fedora 15, AMD Phenom II X4:

sage: f(x) = sin(x)^2/x^2
sage: f.integral(x, -infinity, infinity)
x |--> pi
sage: f.integral(x, -infinity, infinity).n()
3.14159265358979
sage: f.integral(x, -50000, 0).n() + f.integral(x, 0, 50000).n()
3.1415769097886317

everthing fine so far, but

sage: f.integral(x, -50000, 50000).n()
0.06953294323974919
sage: f.integral(x, -5000000, 5000000).n()
0.002749874456609362

Am I over-worked, or is this a bug!?

Georg
Georg

[Andrew Fleckenstein]

I reproduced this in sage 4.8, but with different results:

sage: f(x) = sin(x)^2/x^2
sage: f.integral(x, -50000, 50000).n()
-0.0000200000071537570
sage: f.integral(x, -5000000, 5000000).n()
-2.00000008410959e-7
sage: f.integral(x, -500000000, 500000000).n()
-2.00000000109169e-9

The answer seems to get smaller and smaller as the bounds get farther
apart.

If the bounds get very far apart, something interesting happens:
sage: f.integral(x, -1000000000000000,
1000000000000000)
x |--> -I*gamma(-1, -2000000000000000*I) + I*gamma(-1,
2000000000000000*I) - 1/1000000000000000

I think this has to do with the way sage evaluates an integral. If you
interrupt while it's evaluating, this appears

sage: f.integral(x, -500000001, 500000000).n()
^C---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call
last)

/home/andrew/sage-4.8/<ipython console> in <module>()

/home/andrew/sage-4.8/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._numerical_approx (sage/symbolic/expression.cpp:18004)()

/home/andrew/sage-4.8/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._convert (sage/symbolic/expression.cpp:5089)()

/home/andrew/sage-4.8/local/lib/python2.6/site-packages/sage/functions/other.pyc in _evalf_(self, x, y, parent)
719 """
720 try:
--> 721 return x.gamma_inc(y)
722 except AttributeError:
723 if not (is_ComplexNumber(x)):

/home/andrew/sage-4.8/local/lib/python2.6/site-packages/sage/rings/complex_number.so in sage.rings.complex_number.ComplexNumber.gamma_inc (sage/rings/complex_number.c:12096)()

/home/andrew/sage-4.8/local/lib/python2.6/site-packages/sage/libs/pari/gen.so in sage.libs.pari.gen.gen.incgam (sage/libs/pari/gen.c:19395)()

I have no idea what this means, so maybe someone with more experience
could more accurately diagnose the results.

Cheers!

[Andrew Fleckenstein]

I reproduced this in sage 4.8, but with different results:

sage: f(x) = sin(x)^2/x^2

Cheers!

[kcrisman]

The 4.8 stuff is probably due to a previous Maxima.

In the 5.0 betas, I suspect this particular one is related to things
like

http://trac.sagemath.org/sage_trac/ticket/8321

GSL is hiccuping, I guess.

See also

sage: f.nintegral(x,-50000,50000)

ValueError: Maxima (via quadpack) cannot compute the integral

I don't know if there is any obvious answer to this, based on the
discussion at #8321.

- kcrisman

[Nils Bruin]

On Apr 17, 6:26 pm, kcrisman <kcris...@gmail.com> wrote:
> The 4.8 stuff is probably due to a previous Maxima.
>
> In the 5.0 betas, I suspect this particular one is related to things
> like
>
> http://trac.sagemath.org/sage_trac/ticket/8321

Some of the comments on that ticket also point out that in sage there
doesn't seem to be a way to do arbitrary precision numerical
integration. PARI/GP has some excellent functionality for that and
PARI's intnum accepts its integrand as a black box function, just as
GSL etc.! So with a bit of wrapping, we could use that to provide
arbitrary precision numerical integration too (rather than quadpack/
GSL/numpy/scipy/etc approaches that are all limited to float/double).

Of course, none of the normal gotchas for numerical work magically
disappear, but at least having some hooks to be able to use PARI's
excellent integration routines would be useful to have if more than 53
bits is required. Writing an interface would be a nice project of
fairly limited scope and would require reading up on PARI/GP's
integration routines in order to do it intelligently.