Integration issue

[ketchers]

Sage returns negative value for the integral of a positive function x*cos(x^3) on (0,0.5), if I use abs(cos(x^3))*x, then it gets it correct?

[kcrisman]

This is now http://trac.sagemath.org/sage_trac/ticket/12947.  We've had some issues with incomplete gamma functions translating properly in the past, and/or errors in Maxima, but I didn't have time to either look into that or whether there was another ticket open for this, apologies if there is one - just wanted to make sure this was opened.

- kcrisman

[Robert Dodier]

On 2012-05-14, kcrisman <kcri...@gmail.com> wrote:

> This is now http://trac.sagemath.org/sage_trac/ticket/12947. We've had
> some issues with incomplete gamma functions translating properly in the
> past, and/or errors in Maxima, but I didn't have time to either look into
> that or whether there was another ticket open for this, apologies if there
> is one - just wanted to make sure this was opened.

My first guess is that there is branch cut strangeness going on. Sorry,
I don't have any details. But if you want to investigate, try
integrate(x*cos(x^3), x, 0, u) and then differentiate w.r.t. u, as a
point of departure.

best,

Robert Dodier

[John H Palmieri]

On Sunday, May 13, 2012 8:46:41 PM UTC-7, ketchers wrote:

Sage returns negative value for the integral of a positive function x*cos(x^3) on (0,0.5), if I use abs(cos(x^3))*x, then it gets it correct?

This works for me:

    sage: numerical_integral(x*cos(x^3), 0, 0.5)
    (0.1247560409610376, 1.3850702913602309e-15)

--
John

[JamesHDavenport]

It may be "branch cut strangeness", but if so it is very strange. The integrand is clearly well-behaved, and the integral,
while in terms of the incomplete gamma function, seems to be off the usual branch cut (negative real axis).

[Robert Dodier]

On 2012-05-14, JamesHDavenport <J.H.Da...@bath.ac.uk> wrote:

> It may be "branch cut strangeness", but if so it is very strange. The
> integrand is clearly well-behaved, and the integral,
> while in terms of the incomplete gamma function, seems to be off the usual
> branch cut (negative real axis).

Try domain:complex before calling integrate; that changes the result to
what I think is expected.

I guess (emphasis on guess) that the problem originates not from
gamma_incomplete itself but from terms of the form (-1)^(1/n) which are
the result of simplifying or evaluating gamma_incomplete. Sorry I can't
be more helpful.

best,

Robert Dodier

[Keshav Kini]

John H Palmieri <jhpalm...@gmail.com> writes:
> This works for me:
>
> sage: numerical_integral(x*cos(x^3), 0, 0.5)
> (0.1247560409610376, 1.3850702913602309e-15)

Interesting...

sage: numerical_integral(x*cos(x^3), 0, 0.5)
(0.1247560409610376, 1.3850702913602309e-15)

sage: (x*cos(x^3))(0)
/opt/sage-5.0.rc1/local/lib/python2.7/site-packages/IPython/iplib.py:2260: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...)
exec code_obj in self.user_global_ns, self.user_ns
0

Why does numerical_integral() not trigger the deprecation warning?

-Keshav

----
Join us in #sagemath on irc.freenode.net !

[kcrisman]

> It may be "branch cut strangeness", but if so it is very strange. The
> integrand is clearly well-behaved, and the integral,
> while in terms of the incomplete gamma function, seems to be off the usual
> branch cut (negative real axis).

Try domain:complex before calling integrate; that changes the result to
what I think is expected.

(%i5) display2d:false;

(%o5) false

(%i6) integrate(x*cos(x^3),x);

(%o6) (gamma_incomplete(2/3,%i*x^3)+gamma_incomplete(2/3,-%i*x^3))/6

(%i7) domain:complex;

(%o7) complex

(%i8) integrate(x*cos(x^3),x);

(%o8) ((sqrt(3)*%i-1)*gamma_incomplete(2/3,%i*x^3)

       +(-sqrt(3)*%i-1)*gamma_incomplete(2/3,-%i*x^3))

       *(x^3)^(1/3)

       /(12*x)

But the *definite* integral in both cases is wrong.  Any ideas?

 

(%i1) display2d:false;

(%o1) false

(%i2) integrate(x*cos(x^3),x,0,1/2);

(%o2) gamma_incomplete(2/3,%i/8)/6+gamma_incomplete(2/3,-%i/8)/6-gamma(2/3)/3

(%i3) domain:complex;

(%o3) complex

(%i4) integrate(x*cos(x^3),x,0,1/2);

(%o4) gamma_incomplete(2/3,%i/8)/6+gamma_incomplete(2/3,-%i/8)/6-gamma(2/3)/3

I guess (emphasis on guess) that the problem originates not from
gamma_incomplete itself but from terms of the form (-1)^(1/n) which are
the result of simplifying or evaluating gamma_incomplete. Sorry I can't
be more helpful.

I don't see any of those up here, though, and the gamma_incomplete evaluation is correct (gives the same via W|A, Sage = Pari in my version, mpmath, and Maxima).  I think that Maxima is somehow using the "real" antiderivative, if that makes sense - is that possible, Robert?

[kcrisman]

The same reason that plot and integral don't, because we're not "calling" them in the same way.  It makes sense to integrate symbolic expressions and to plot them.  It's true that we need to unify our integration command syntax (see http://trac.sagemath.org/sage_trac/ticket/7763).

- kcrisman 

[Robert Dodier]

On 2012-05-15, kcrisman <kcri...@gmail.com> wrote:

> (%i3) domain:complex;
>
> (%o3) complex
> (%i4) integrate(x*cos(x^3),x,0,1/2);
>
> (%o4)
> gamma_incomplete(2/3,%i/8)/6+gamma_incomplete(2/3,-%i/8)/6-gamma(2/3)/3

Hmm. I get a different result. I am using the current Git version.

domain : complex;
integrate (x*cos(x^3), x, 0, 1/2);
=>
%i*gamma_incomplete(2/3,%i/8)/(4*sqrt(3))
-gamma_incomplete(2/3,%i/8)/12-%i*gamma_incomplete(2/3,-%i/8)/(4*sqrt(3))
-gamma_incomplete(2/3,-%i/8)/12+gamma(2/3)/6
expand (float (%));
=> .1247560409610377

That's gratifying, but the problem, as I'm sure you know, is that the
user won't know they have to change a global variable.

> I don't see any of those up here, though, and the gamma_incomplete
> evaluation is correct (gives the same via W|A, Sage = Pari in my version,
> mpmath, and Maxima). I think that Maxima is somehow using the "real"
> antiderivative, if that makes sense - is that possible, Robert?

It seems plausible, but I don't know the integration code very well.

best

Robert Dodier

[kcrisman]

> (%i3) domain:complex;
>
> (%o3) complex
> (%i4) integrate(x*cos(x^3),x,0,1/2);
>
> (%o4)
> gamma_incomplete(2/3,%i/8)/6+gamma_incomplete(2/3,-%i/8)/6-gamma(2/3)/3

Hmm. I get a different result. I am using the current Git version.

Great, I didn't realize some code had changed - I get this same thing below in 5.27.0.

 

domain : complex;
integrate (x*cos(x^3), x, 0, 1/2);
 =>
%i*gamma_incomplete(2/3,%i/8)/(4*sqrt(3))
 -gamma_incomplete(2/3,%i/8)/12-%i*gamma_incomplete(2/3,-%i/8)/(4*sqrt(3))
 -gamma_incomplete(2/3,-%i/8)/12+gamma(2/3)/6
expand (float (%));
 => .1247560409610377

That's gratifying, but the problem, as I'm sure you know, is that the
user won't know they have to change a global variable.

If all integrals still work with domain:complex, we could conceivably set this in the integration code.  However, we *already* set `domain:complex` in the startup code for the "calculus" copy of Maxima for exactly this reason... so I guess that this would be fixed by upgrading Maxima?

[Keshav Kini]

kcrisman <kcri...@gmail.com> writes:
> On Tuesday, May 15, 2012 2:36:34 AM UTC-4, Keshav Kini wrote:
> Why does numerical_integral() not trigger the deprecation warning?
>
> The same reason that plot and integral don't, because we're not "calling" them
> in the same way. It makes sense to integrate symbolic expressions and to plot
> them.

I don't know if I agree. I mean, certainly it makes sense to integrate
symbolic expressions and to plot them, but the point of the deprecation
warning is that you should specify which variable you want to integrate
with respect to, or plot with respect to.

Then again, with plot(), if there's more than one free variable in the
symbolic expression, it doesn't make sense to plot it anyway. The same
is true of numerical integration (though not of symbolic integration).
So I guess I answered my own question.

Note: the above reasoning doesn't apply to plot3d() - there is still
ambiguity as to which of two free variables should be on the x axis and
which on the y axis. And maybe that's why plot3d(), unlike plot(), does
seem to generate the deprecation warning.

[Jason Grout]

On 5/15/12 8:33 PM, Keshav Kini wrote:
> And maybe that's why plot3d(), unlike plot(), does
> seem to generate the deprecation warning.

Sorry---what plot command doesn't generate a deprecation warning?

Thanks,

Jason

[Keshav Kini]

For example:

sage: plot(x^2, (0, 1))

sage:

[Jason Grout]

On 5/15/12 9:17 PM, Keshav Kini wrote:
> plot(x^2, (0, 1))

I definitely think that should give a deprecation warning (I think I've
been advocating for that to give a deprecation warning for a long time).
For example, I think this is confusing:

plot(x^2+y-x^2,(0,1))

Jason

[kcrisman]

Yes yes, we've had this conversation before.  Then maybe we shouldn't simplify such things automatically, in which case a warning would actually make sense - who actually types in "x^2+y-x^2" to plot it in one variable?  Like Keshav said, this answers its own question - I can't remember ever having had an actual support request on this issue.   On the most basic commands, requiring extra syntactical sugar is a recipe for stultification, and plot is surely one of those commands.  I really don't think we need to get that language legalistic on things like this.

[Jason Grout]

The point is that I can't tell if a function f is actually a function of
one variable because simplification is Hard. The example above is
contrived, but illustrates the point that we may have no idea what is
being plotted. Explicit is better than implicit, etc.

You're right that we've had I don't know how many emails and threads
over this issue. I guess I still haven't changed my mind. One big
thing for making the syntax above illegal is that right now there is a
very ugly wart in the code special-casing that case, and it is
inconsistent with how other things work (for example, note the confusion
on this thread).

Anyway, I'm not going to do anything about it in the near future, so
I'll note my wish that the above had a deprecation warning and go back
to grading.

Thanks,

Jason

[Jason Grout]

On 5/15/12 9:42 PM, Jason Grout wrote:

> Anyway, I'm not going to do anything about it in the near future, so
> I'll note my wish that the above had a deprecation warning and go back
> to grading.

At the very least, we should maybe print out a message saying that we
are going to assume that the horizontal axis represents the variable x.

Thanks,

Jason

[ketchers]

I don't know how to get sage to understand "domain : complex" so I tried with assume and here is what happened. Does it make sense? 

[kcrisman]

On Thursday, May 17, 2012 1:52:38 AM UTC-4, ketchers wrote:

I don't know how to get sage to understand "domain : complex" so I tried with assume and here is what happened. Does it make sense? 

Yes, it does.  Our assumptions go through Maxima, and apparently assuming a variable is complex does the job.

sage: a = integrate(x*cos(x^3),(x,0,0.5)).n()

sage: assume(x,'complex')

sage: b = integrate(x*cos(x^3),(x,0,0.5)).n()

sage: assumptions()

[x is complex]

But apparently we didn't think of this redundancy.  I suppose it's possible for a variable to be real and complex, but probably this is not what is intended, as your post implies.   (Try to actually cut and paste commands, having it in the image makes it tedious to reproduce.)

sage: assume(x,'real')

sage: assumptions()

[x is complex, x is real]

sage: maxima_calculus('domain')

complex

sage: forget()

sage: c = integrate(x*cos(x^3),(x,0,0.5)).n()

sage: a; b; c

-0.0677842754623305

0.1247560409610376

-0.0677842754623305

 

[kcrisman]

Yes, it does.  Our assumptions go through Maxima, and apparently assuming a variable is complex does the job.

sage: a = integrate(x*cos(x^3),(x,0,0.5)).n()

sage: assume(x,'complex')

sage: b = integrate(x*cos(x^3),(x,0,0.5)).n()

sage: assumptions()

[x is complex]

Actually,

sage: assume(x,'complex')

sage: integrate(x*cos(x^3),(x,0,0.5))

integrate(x*cos(x^3), x, 0, 0.5)

so the reason it works is because Maxima doesn't even evaluate it in the version of Maxima we have, so it goes straight to numerical integration with .n(). 

Anyway, upgrading our Maxima will take care of this integral.