Strange behavior when integrating 1/x^p

[Aaron Tresham]

I'm assuming the problem I'm encountering is not in SMC, but in whatever program is being used to integrate. But perhaps someone here knows why this happens.

If I try the following integrals, the first three work, but the last returns an error (RuntimeError: ECL says: CEXPT only defined for non-negative integral exponents.).

integral(1/x^(1/2),x,1,10)

integral(1/x^.5,x,1,10)

integral(1/x^(3/2),x,1,10)

integral(1/x^1.5,x,1,10)

If the power is pi or sqrt(2), it also works fine. And if I change all the powers to negative, then it works.

So it seems non-integer rational numbers greater than 1 must be written as a fraction.

Thanks,

Aaron

[Harald Schilly]

This is a problem within Maxima.
I am not sure what you are trying to do, there is also "integrate"
(not "integral") and "numerical_integral". Besides that, you can
change the back-end, which is used to integrate, to sympy:

sage: integrate(1/x^1.5,x,1,10, algorithm="sympy")
1.36754446796632

-- harald

> --
> You received this message because you are subscribed to the Google Groups
> "sage-cloud" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sage-cloud+...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sage-cloud/700ac3b2-1e3a-4e34-8679-6e38258ac7fa%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[William Stein]

(Response from the lead Maxima developer, where Maxima is what is
called by Sage to do this computation.)


---------- Forwarded message ----------
From: Robert Dodier <robert...@gmail.com>
Date: Wed, Oct 15, 2014 at 2:19 PM
Subject: [sage-support] Re: Fwd: [sage-cloud] Strange behavior when
integrating 1/x^p
To: sage-s...@googlegroups.com

On 2014-10-15, William Stein <wst...@gmail.com> wrote:

> If I try the following integrals, the first three work, but the last
> returns an error (RuntimeError: ECL says: CEXPT only defined for
> non-negative integral exponents.).
>
> integral(1/x^(1/2),x,1,10)
> integral(1/x^.5,x,1,10)
> integral(1/x^(3/2),x,1,10)
> integral(1/x^1.5,x,1,10)

The immediate cause of the error is that Sage enables the keepfloat
flag, otherwise Maxima converts floats to rationals when trying to
compute integrals. Some (maybe a lot) of the code for working with
polynomials and rational functions expects only integers and rationals,
so the presence of floats causes trouble.

I suppose that we have to decide between fixing up Maxima to handle
floats everywhere that keepfloat might introduce them, or do away
with keepfloat. I'm in favor of the latter.

Probably would be helpful if someone would file a bug report.
http://sourceforge.net/p/maxima/bugs

best,

Robert Dodier

--
You received this message because you are subscribed to the Google

Groups "sage-support" group.

To unsubscribe from this group and stop receiving emails from it, send

an email to sage-support...@googlegroups.com.
To post to this group, send email to sage-s...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.

For more options, visit https://groups.google.com/d/optout.

--
William Stein
Professor of Mathematics
University of Washington
http://wstein.org

[Aaron Tresham]

Thank you, Harald. Using sympy solves the problem for me.

I did notice some strange behavior with the improper integrals below when using sympy. The first returns +Infinity, but the second returns an error (CoeffExpValueError: expr not of form a*x**b)

integral(1/x^(1/2),x,1,+Infinity,algorithm='sympy')

integral(1/x^(.5),x,1,+Infinity,algorithm='sympy')

I don't encounter this difference if I replace the +Infinity with a number, or if the power is bigger than 1 (so the integral converges).

If I use Maxima, then both return ValueError: Integral is divergent.

[Harald Schilly]

On Thu, Oct 16, 2014 at 2:05 AM, Aaron Tresham <atre...@gmail.com> wrote:
> (CoeffExpValueError: expr not of form a*x**b)
>

> integral(1/x^(.5),x,1,+Infinity,algorithm='sympy')

I think the core problem is, that .5 is inherently inexact -- but --
you are asking for exact solutions. So, that behavior is actually not
bad, at least from my POV.
Alternatively, there are various quadrature routines for inexact
integration. They are e.g. in mpmath, which is in Sage and in Sympy.
There are also quadrature formulas in Sage, but well,
numerical_integral also doesn't do what I want. Maybe there is another
one.

Anyways, as you can see, this doesn't work well for pure numerical integration:

import mpmath as mp

>>> import mpmath as mp

>>> mp.quad(lambda x : 1/x^(.5), [1,+Infinity])
mpf('7464407990.9940443')

>>> mp.quad(lambda x : 1/x^(.5), [1,+Infinity], maxdegree=5)
mpf('6190466468.2364979')

>>> mp.quad(lambda x : 1/x^(.5), [1,+Infinity], maxdegree=10)
mpf('8636927446.1778164')