Bug in integrating x*cos(x^3) #12947

kcrisman · 2012-05-14T12:33:57Z

sage: numerical_integral(x*cos(x^3),0,.5)
(0.1247560409610376, 1.3850702913602309e-15)
sage: integrate(x*cos(x^3),(x,0,1/2)).n()
-0.0677842754623305

Given

sage: integrate(x*cos(x^3),x)
1/12*((-I*sqrt(3) - 1)*gamma(2/3, -I*x^3) + (I*sqrt(3) - 1)*gamma(2/3, I*x^3))*(x^3)^(1/3)/x
sage: integrate(x*cos(x^3),(x,0,1/2))
-1/3*gamma(2/3) + 1/6*gamma(2/3, -1/8*I) + 1/6*gamma(2/3, 1/8*I)

this is probably something where Maxima is returning too many imaginary things and something isn't cancelling right either there or in Sage proper? We've had trouble with either side of this with incomplete gamma functions before.

CC: @orlitzky

Component: calculus

Author: Peter Bruin

Branch/Commit: 72706eb

Reviewer: Ralf Stephan

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

The text was updated successfully, but these errors were encountered:

kcrisman · 2012-05-15T02:30:47Z

comment:1

sage: assume(u>0)
sage: integrate(x*cos(x^3),x,0,u)
-1/3*gamma(2/3) + 1/6*gamma(2/3, -I*u^3) + 1/6*gamma(2/3, I*u^3)

To compare the "right" answer, do

sage: f(u) = integrate(x*cos(x^3),x,0,u)
sage: plot(f,(u,0,2*pi))
sage: plot(lambda u: numerical_integral(x*cos(x^3),0,u)[0],(u,0,2*pi))

Note how they don't really look quite similar, not really a branch cut thing? Also, in Maxima

(%i2) display2d:false;
(%o2) false
(%i3) f:integrate(x*cos(x^3),x);
(%o3) (gamma_incomplete(2/3,%i*x^3)+gamma_incomplete(2/3,-%i*x^3))/6
(%i5) diff(f,x);
(%o5) (-3*%i*x*%e^(%i*x^3)/(-1)^(1/6)-3*%i*x*%e^-(%i*x^3)/(-1)^(1/6))/6

kcrisman · 2012-05-15T14:04:13Z

comment:2

Robert Dodier of Maxima suggests the following later in the thread.


(%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 this doesn't resolve the original issue.


(%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

kcrisman · 2012-05-15T15:46:01Z

comment:3

Robert points out that in 5.27.0 this is at least a little better, though we still need to change domain.


Maxima 5.27.0 http://maxima.sourceforge.net
using Lisp SBCL 1.0.24
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) display2d:false;

(%o1) false
(%i2) domain:complex;

(%o2) complex
(%i3) integrate (x*cos(x^3), x, 0, 1/2); 

(%o3) %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
(%i4) expand(float(%));

(%o4) .1247560409610377

kcrisman · 2012-05-15T15:50:12Z

comment:4

But since we already use domain:complex in the initialization code for Maxima for integrals and such, presumably an upgrade would fix this?

kcrisman · 2012-05-15T15:50:12Z

Upstream: Fixed upstream, in a later stable release.

orlitzky · 2012-05-23T17:40:09Z

comment:5

See my last comment on #11238 for why we can't have nice things.

That's the only real issue with maxima-5.27. Some messages are getting mangled e.g. when maxima asks us for more information, but they aren't substantial regressions.

pjbruin · 2014-05-30T22:49:36Z

comment:9

This seems to have been fixed some time ago; it works correctly in both Sage 5.13 and 6.2.

sage: integrate(x*cos(x^3),(x,0,1/2)).n()
0.124756040961038

pjbruin · 2014-05-30T22:49:36Z

Changed upstream from Fixed upstream, in a later stable release. to none

pjbruin · 2014-05-30T22:49:36Z

Work Issues: add doctest

pjbruin · 2014-05-30T23:09:20Z

Changed work issues from add doctest to none

pjbruin · 2014-05-30T23:09:20Z

Branch: u/pbruin/12947-maxima_integral

pjbruin · 2014-05-30T23:09:20Z

Author: Peter Bruin

pjbruin · 2014-05-30T23:09:20Z

Commit: 35f6081

rwst · 2014-05-31T08:13:27Z

comment:11

Funny I get:

File "src/sage/interfaces/maxima_lib.py", line 736, in sage.interfaces.maxima_lib.MaximaLib.sr_integral
Failed example:
    integrate(x*cos(x^3),(x,0,1/2)).n()
Expected:
    0.124756040961038
Got:
    0.124756040961038 + 2.77555756156289e-17*I

Please set positive when you think this is resolved.

rwst · 2014-05-31T08:13:27Z

Reviewer: Ralf Stephan

sagetrac-git · 2014-05-31T10:04:39Z

Branch pushed to git repo; I updated commit sha1. New commits:

`3d9a07c`	`Trac 12947: add tolerance to doctest`

sagetrac-git · 2014-05-31T10:04:39Z

Changed commit from 35f6081 to 3d9a07c

pjbruin · 2014-05-31T10:06:08Z

comment:13

On the system I tested this on (GNU/Linux x86_64) this gives the expected answer, both with and without the Maxima upgrade of #13973. I just checked on ARM; the first time I got the same answer as you, but then I merged this branch (did not even rebuild) and only got the expected answer from then on. I added # abs tol 1e-16; can you check if this solves the problem for you?

kcrisman · 2014-05-31T13:02:26Z

comment:14

I would argue that the complex bit means that we are not quite ready to merge this if it happens on some platforms. The question is whether this is a Maxima bug or a Sage/Pynac simplification bug.

rwst · 2014-05-31T13:08:00Z

comment:15

Would you equally mind if there is a real digit wrong in the 17-th place? Regardless, it's fine on AMD now.

rwst · 2014-05-31T13:11:24Z

comment:16

There are 60(!) lines in the sage code with # abs tol and nearly all have lower tolerance (in the sense of 'larger error')

pjbruin · 2014-05-31T13:27:33Z

comment:17

The complex bit apparently comes from applying Pynac's evalf() to the following expression:

sage: e = integrate(x*cos(x^3),(x,0,1/2)); e
1/12*I*sqrt(3)*gamma(2/3, 1/8*I) - 1/12*I*sqrt(3)*gamma(2/3, -1/8*I) - 1/12*gamma(2/3, 1/8*I) - 1/12*gamma(2/3, -1/8*I) + 1/6*gamma(2/3)
sage: e.n(53)
0.124756040961038 + 2.77555756156289e-17*I
sage: e._convert({'parent': CC})  # called by .n()
0.124756040961038 + 2.77555756156289e-17*I

(This is on ARM; the imaginary part does not appear on x86_64.)

There are known precision problems in the incomplete Gamma function, so this may be related to #7099, although we only want the default 53 bits of precision here.

pjbruin · 2014-05-31T15:48:20Z

comment:18

The presence of the imaginary part does not depend on #7099. At least on ARM it is essentially random (maybe it depends on FPU rounding flags or something similar):

$ ./sage -c 'print integrate(x*cos(x^3),(x,0,1/2)).n()'
0.124756040961038
$ ./sage -c 'print integrate(x*cos(x^3),(x,0,1/2)).n()'
0.124756040961038 - 2.77555756156289e-17*I
$ ./sage -c 'print integrate(x*cos(x^3),(x,0,1/2)).n()'
0.124756040961038 - 1.38777878078145e-17*I
$ ./sage -c 'print integrate(x*cos(x^3),(x,0,1/2)).n()'
0.124756040961038 + 1.38777878078145e-17*I

kcrisman · 2014-06-01T01:48:35Z

comment:19

Would you equally mind if there is a real digit wrong in the 17-th place? Regardless, it's fine on AMD now.

Haha! No, but as I said sometimes the complex bit indicates a different kind of problem, or has in the past. Given Peter's experiment in comment:18, though, I guess we can make that a different ticket.

pjbruin · 2014-06-02T00:12:44Z

comment:20

The imaginary part turns out to be "simply" the fact that floating-point addition is not associative, and for some reason the terms in the expression e from comment:17 are added in random order by .n().

sage: e = integrate(x*cos(x^3),(x,0,1/2)); e
1/12*I*sqrt(3)*gamma(2/3, 1/8*I) - 1/12*I*sqrt(3)*gamma(2/3, -1/8*I) - 1/12*gamma(2/3, 1/8*I) - 1/12*gamma(2/3, -1/8*I) + 1/6*gamma(2/3)
sage: v = [term.n() for term in e.operands()]; v
[0.0454319516149362 + 0.166098636946833*I,
 0.0454319516149362 - 0.166098636946833*I,
 -0.0958970927532841 + 0.0262301494946934*I,
 -0.0958970927532841 - 0.0262301494946934*I,
 0.225686323237733]
sage: uniq([sum(p.action(v)) for p in permutations(5)])
[0.124756040961038 - 2.77555756156289e-17*I,
 0.124756040961038,
 0.124756040961038 + 2.77555756156289e-17*I,
 0.124756040961038 - 1.38777878078145e-17*I,
 0.124756040961038,
 0.124756040961038 + 1.38777878078145e-17*I,
 0.124756040961038 - 2.77555756156289e-17*I,
 0.124756040961038,
 0.124756040961038 + 2.77555756156289e-17*I,
 0.124756040961038 - 1.38777878078145e-17*I,
 0.124756040961038 + 1.38777878078145e-17*I,
 0.124756040961038 - 2.77555756156289e-17*I,
 0.124756040961038 + 2.77555756156289e-17*I,
 0.124756040961038]

This gives 14 different outcomes depending on the summation order (some of the differences are hidden due to the number of digits printed). I don't think this can be solved easily, so I propose we stick to the # abs tol fix.

sagetrac-git · 2014-06-11T22:37:53Z

Branch pushed to git repo; I updated commit sha1. New commits:

`8274e47`	`Merge branch 'develop' into ticket/12947-numerical_integral`

sagetrac-git · 2014-06-11T22:37:53Z

Changed commit from 3d9a07c to 8274e47

rwst · 2014-06-29T08:59:44Z

comment:22

As already said ...

vbraun · 2014-06-30T13:12:23Z

comment:23

The # abs tol just changes how the floating point numbers are compared, it doesn't know about complex numbers.

sage -t --long src/sage/interfaces/maxima_lib.py
**********************************************************************
File "src/sage/interfaces/maxima_lib.py", line 763, in sage.interfaces.maxima_lib.MaximaLib.sr_integral
Failed example:
    integrate(x*cos(x^3),(x,0,1/2)).n()  # abs tol 1e-16
Expected:
    0.124756040961038
Got:
    0.124756040961038 - 1.38777878078145e-17*I
**********************************************************************

and on a different machine

sage -t --long src/sage/interfaces/maxima_lib.py
**********************************************************************
File "src/sage/interfaces/maxima_lib.py", line 763, in sage.interfaces.maxima_lib.MaximaLib.sr_integral
Failed example:
    integrate(x*cos(x^3),(x,0,1/2)).n()  # abs tol 1e-16
Expected:
    0.124756040961038
Got:
    0.124756040961038 - 2.77555756156289e-17*I
**********************************************************************

sagetrac-git · 2014-07-18T15:25:44Z

Changed commit from 8274e47 to 72706eb

sagetrac-git · 2014-07-18T15:25:44Z

Branch pushed to git repo; I updated commit sha1. New commits:

`72706eb`	`Trac 12947: handle numerical noise correctly`

pjbruin · 2014-07-18T15:29:35Z

comment:26

This should solve it; probably no reason to ask for another review.

vbraun · 2014-07-19T04:57:53Z

Changed branch from u/pbruin/12947-maxima_integral to 72706eb

kcrisman added this to the sage-5.11 milestone May 14, 2012

kcrisman added c: calculus labels May 14, 2012

kcrisman assigned burcin May 14, 2012

jdemeyer modified the milestones: sage-5.11, sage-5.12 Aug 13, 2013

sagetrac-vbraun-spam mannequin modified the milestones: sage-6.1, sage-6.2 Jan 30, 2014

sagetrac-vbraun-spam mannequin modified the milestones: sage-6.2, sage-6.3 May 6, 2014

pjbruin added p: trivial / 5 and removed p: major / 3 labels May 30, 2014

pjbruin added the s: needs review label May 30, 2014

rwst added s: positive review and removed s: needs review labels Jun 29, 2014

vbraun added s: needs work and removed s: positive review labels Jun 30, 2014

pjbruin added s: positive review and removed s: needs work labels Jul 18, 2014

vbraun removed the s: positive review label Jul 19, 2014

vbraun closed this as completed in dc8130a Jul 19, 2014

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Bug in integrating x*cos(x^3) #12947

Bug in integrating x*cos(x^3) #12947

kcrisman commented May 14, 2012

kcrisman commented May 15, 2012

kcrisman commented May 15, 2012

kcrisman commented May 15, 2012

kcrisman commented May 15, 2012

kcrisman commented May 15, 2012

orlitzky commented May 23, 2012

pjbruin commented May 30, 2014

pjbruin commented May 30, 2014

pjbruin commented May 30, 2014

pjbruin commented May 30, 2014

pjbruin commented May 30, 2014

pjbruin commented May 30, 2014

pjbruin commented May 30, 2014

rwst commented May 31, 2014

rwst commented May 31, 2014

sagetrac-git mannequin commented May 31, 2014

sagetrac-git mannequin commented May 31, 2014

pjbruin commented May 31, 2014

kcrisman commented May 31, 2014

rwst commented May 31, 2014

rwst commented May 31, 2014

pjbruin commented May 31, 2014

pjbruin commented May 31, 2014

kcrisman commented Jun 1, 2014

pjbruin commented Jun 2, 2014

sagetrac-git mannequin commented Jun 11, 2014

sagetrac-git mannequin commented Jun 11, 2014

rwst commented Jun 29, 2014

vbraun commented Jun 30, 2014

sagetrac-git mannequin commented Jul 18, 2014

sagetrac-git mannequin commented Jul 18, 2014

pjbruin commented Jul 18, 2014

vbraun commented Jul 19, 2014

Bug in integrating x*cos(x^3) #12947

Bug in integrating x*cos(x^3) #12947

Comments

kcrisman commented May 14, 2012

kcrisman commented May 15, 2012

kcrisman commented May 15, 2012

kcrisman commented May 15, 2012

kcrisman commented May 15, 2012

kcrisman commented May 15, 2012

orlitzky commented May 23, 2012

pjbruin commented May 30, 2014

pjbruin commented May 30, 2014

pjbruin commented May 30, 2014

pjbruin commented May 30, 2014

pjbruin commented May 30, 2014

pjbruin commented May 30, 2014

pjbruin commented May 30, 2014

rwst commented May 31, 2014

rwst commented May 31, 2014

sagetrac-git mannequin commented May 31, 2014

sagetrac-git mannequin commented May 31, 2014

pjbruin commented May 31, 2014

kcrisman commented May 31, 2014

rwst commented May 31, 2014

rwst commented May 31, 2014

pjbruin commented May 31, 2014

pjbruin commented May 31, 2014

kcrisman commented Jun 1, 2014

pjbruin commented Jun 2, 2014

sagetrac-git mannequin commented Jun 11, 2014

sagetrac-git mannequin commented Jun 11, 2014

rwst commented Jun 29, 2014

vbraun commented Jun 30, 2014

sagetrac-git mannequin commented Jul 18, 2014

sagetrac-git mannequin commented Jul 18, 2014

pjbruin commented Jul 18, 2014

vbraun commented Jul 19, 2014