Have numerical evaluation of unevaluated integrals call numerical integral

[dandrake]

----------------------------------------------------------------------
| SAGE Version 3.0.6, Release Date: 2008-07-30                       |
| Type notebook() for the GUI, and license() for information.        |
----------------------------------------------------------------------

sage: x = var('x')
sage: integrate(x^2.7 * e^(-2.4*x), x, 0, 3).n()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

/home/drake/sage-3.0.6.final/<ipython console> in <module>()

/home/was/s/local/lib/python2.5/site-packages/sage/calculus/calculus.py in numerical_approx(self, prec, digits)
   1266         except TypeError:
   1267             # try to return a complex result
-> 1268             approx = self._complex_mpfr_field_(ComplexField(prec))
   1269 
   1270         return approx

/home/was/s/local/lib/python2.5/site-packages/sage/calculus/calculus.py in _complex_mpfr_field_(self, field)
   1419 
   1420     def _complex_mpfr_field_(self, field):
-> 1421         raise TypeError
   1422 
   1423     def _complex_double_(self, C):

TypeError: 
sage: 

Oddly, the plot function has no difficulty, so some part of Sage can numerically evaluate the function:

plot(x^2.7 * e^(-2.4*x), x, 0, 3)

works fine.

Some values for the exponents do work -- it seems like the exponent of x needs to be an integer or half-integer:

(2.7, -2.4): this is the above example
(27/10, -2.4): same error as above
(1.5, -2.4): works
(1.6, -2.4): same error as above
(1.6, -2.0): same error as above
(1.0, -2.4): works
(5.5, -2.4): works

Component: calculus

Keywords: integration integral calculus symbolic numerical

Author: Golam Mortuza Hossain, Burcin Erocal

Reviewer: Dan Drake

Merged: sage-4.4.alpha0

Issue created by migration from https://trac.sagemath.org/ticket/3863

[dandrake]

comment:1

With the new symbolics in 4.0.rc0, this still doesn't work; the error is now:

sage: integrate(x^2.7 * exp(-2.4*x), x, 0, 3).n()
ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (10003, 0))

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

/home/drake/.sage/temp/sage.math.washington.edu/22166/_home_drake__sage_init_sage_0.py in <module>()

/scratch/drake/4.0.rc0/local/lib/python2.5/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.n (sage/symbolic/expression.cpp:15691)()

/scratch/drake/4.0.rc0/local/lib/python2.5/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.pyobject (sage/symbolic/expression.cpp:2641)()

TypeError: self must be a numeric expression

Putting parentheses around the 2.7 and -2.4 didn't change anything.

[kcrisman]

comment:2

The essential problem is that with half-integer or integer exponents of x, Maxima can symbolically integrate this using Erf. Otherwise it can't (that doesn't mean it's not possible, just that Maxima doesn't know). Before, Sage tried to turn the expression into a complex one if it couldn't evaluate it, but that doesn't do much for a real (unevaluated) integral; now the new symbolics just complain that it's not numeric, which of course it isn't.

My view is that the correct fix is to put some kind of check in for when "integrate" is part of the output string fed into .n() and in that case at least attempt to use numerical_integral or something. Of course that has the problem that things like

sage: integrate(1/(1+x^7))
1/7*log(x + 1) - 1/7*integrate((x^5 - 2*x^4 + 3*x^3 - 4*x^2 + 5*x - 6)/(x^6 - x^5 + x^4 - x^3 + x^2 - x + 1), x)

Maxima is able to partly symbolically integrate, leaving

sage: integrate(1/(1+x^7),x,0,1)

in trouble. But at least a check for "integrate" at the beginning could help.

[dandrake]

comment:3

Replying to @kcrisman:

The essential problem is that with half-integer or integer exponents of x, Maxima can symbolically integrate this using Erf. Otherwise it can't (that doesn't mean it's not possible, just that Maxima doesn't know). Before, Sage tried to turn the expression into a complex one if it couldn't evaluate it, but that doesn't do much for a real (unevaluated) integral; now the new symbolics just complain that it's not numeric, which of course it isn't.

Ah, that makes sense. I don't mind that just naively running .n() doesn't work; we could have, like Mathematica, some sort of numerical_integrate that would try to evaluate the integral, and then punt to something like Simpson's rule to just estimate it. I'm sure, though, that there's vastly better ways to actually do numerical integrals than just a plain vanilla Calculus 1 application of Simpson's Rule.

[kcrisman]

comment:4

Oh, totally, and of course we have numerical_integrate and nintegral (one is Maxima, one is GSL I think). So I guess what I'm saying is that when someone has time (which won't be me until June due to the (US) holiday and then a conference) it should be relatively straightforward to at least do something in the generic code for .n() on whatever kind of object (symbolic?) an expression like integrate(f(x),x,a,b) is like

try:
    usual code
except TypeError:
    if "integrate" is in the string:
        either reraise the TypeError with a message suggesting the use of numerical_integrate
        or actually try to string-magic replace "integrate" with "numerical_integrate", which should be possible (you'd have to take [0] of that result, of course)
    else:
        do whatever used to happen with a TypeError

I think that it is reasonable for a user to expect that .n() numerically evaluates an expression as best as possible, even one Maxima can't evaluate! And all these methods are much more sophisticated than Simpson, certainly, so that's not the issue.

[mwhansen]

comment:5

Actually, the solution is much cleaner since .n() will call the _evalf_ method on the integrate function which we can just have resort to a numerical integral if it is definite and has no free variables.

I'll look into this later.

[kcrisman]

Changed keywords from integration integral calculus symbolic to integration integral calculus symbolic numerical

[kcrisman]

comment:6

Update: in 4.3.alpha1 with Maxima 5.20.1, we now get

sage: integrate(x^2.7 * e^(-2.4*x), x, 0, 3)
119/6144*2^(3/5)*3^(3/10)*5^(7/10)*gamma(7/10) - 125/20736*2^(3/5)*3^(3/10)*5^(7/10)*gamma_incomplete(37/10, 36/5)

But it won't evaluate the gamma_incomplete, since for some reason we aren't translating it to gamma_inc or incomplete_gamma, which are the supported functions; however, otherwise it is correct (as comparing with results of numerical_integral).

This does not fix the problem, of course, but I will change the summary to get at the fundamental thing mhansen and I discussed, and open a separate ticket (if it doesn't exist) for the gamma_incomplete not being translated correctly from Maxima. That is #7748.

Do we in the meantime have the _evalf_ method on a symbolic integral that can be changed?