Opened 10 years ago
Closed 8 years ago
#12947 closed defect (fixed)
Bug in integrating x*cos(x^3)
Reported by:  kcrisman  Owned by:  burcin 

Priority:  trivial  Milestone:  sage6.3 
Component:  calculus  Keywords:  
Cc:  mjo  Merged in:  
Authors:  Peter Bruin  Reviewers:  Ralf Stephan 
Report Upstream:  N/A  Work issues:  
Branch:  72706eb (Commits, GitHub, GitLab)  Commit:  72706eb6697fba361e138ddaca0765f8415eb8f3 
Dependencies:  Stopgaps: 
Description
From this sagesupport thread:
sage: numerical_integral(x*cos(x^3),0,.5) (0.1247560409610376, 1.3850702913602309e15) 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.
Change History (27)
comment:1 Changed 10 years ago by
comment:2 Changed 10 years ago by
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)*%i1)*gamma_incomplete(2/3,%i*x^3) +(sqrt(3)*%i1)*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)/6gamma(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)/6gamma(2/3)/3
comment:3 Changed 10 years ago by
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
comment:4 Changed 10 years ago by
 Cc mjo added
 Report Upstream changed from N/A to Fixed upstream, in a later stable release.
But since we already use domain:complex
in the initialization code for Maxima for integrals and such, presumably an upgrade would fix this?
comment:5 Changed 10 years ago by
See my last comment on #11238 for why we can't have nice things.
That's the only real issue with maxima5.27. Some messages are getting mangled e.g. when maxima asks us for more information, but they aren't substantial regressions.
comment:6 Changed 9 years ago by
 Milestone changed from sage5.11 to sage5.12
comment:7 Changed 8 years ago by
 Milestone changed from sage6.1 to sage6.2
comment:8 Changed 8 years ago by
 Milestone changed from sage6.2 to sage6.3
comment:9 Changed 8 years ago by
 Priority changed from major to trivial
 Report Upstream changed from Fixed upstream, in a later stable release. to N/A
 Work issues set to add doctest
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
comment:10 Changed 8 years ago by
 Branch set to u/pbruin/12947maxima_integral
 Commit set to 35f60813c540f3eaabc45200f3e61f3a227ad1d2
 Status changed from new to needs_review
 Work issues add doctest deleted
comment:11 Changed 8 years ago by
 Reviewers set to Ralf Stephan
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.77555756156289e17*I
Please set positive when you think this is resolved.
comment:12 Changed 8 years ago by
 Commit changed from 35f60813c540f3eaabc45200f3e61f3a227ad1d2 to 3d9a07cd45fd002020a1f67cb2394af5a38d76a0
Branch pushed to git repo; I updated commit sha1. New commits:
3d9a07c  Trac 12947: add tolerance to doctest

comment:13 Changed 8 years ago by
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 1e16
; can you check if this solves the problem for you?
comment:14 Changed 8 years ago by
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.
comment:15 followup: ↓ 19 Changed 8 years ago by
Would you equally mind if there is a real digit wrong in the 17th place? Regardless, it's fine on AMD now.
comment:16 Changed 8 years ago by
There are 60(!) lines in the sage code with # abs tol
and nearly all have lower tolerance (in the sense of 'larger error')
comment:17 Changed 8 years ago by
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.77555756156289e17*I sage: e._convert({'parent': CC}) # called by .n() 0.124756040961038 + 2.77555756156289e17*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.
comment:18 Changed 8 years ago by
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.77555756156289e17*I $ ./sage c 'print integrate(x*cos(x^3),(x,0,1/2)).n()' 0.124756040961038  1.38777878078145e17*I $ ./sage c 'print integrate(x*cos(x^3),(x,0,1/2)).n()' 0.124756040961038 + 1.38777878078145e17*I
comment:19 in reply to: ↑ 15 Changed 8 years ago by
Would you equally mind if there is a real digit wrong in the 17th 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.
comment:20 Changed 8 years ago by
The imaginary part turns out to be "simply" the fact that floatingpoint 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.77555756156289e17*I, 0.124756040961038, 0.124756040961038 + 2.77555756156289e17*I, 0.124756040961038  1.38777878078145e17*I, 0.124756040961038, 0.124756040961038 + 1.38777878078145e17*I, 0.124756040961038  2.77555756156289e17*I, 0.124756040961038, 0.124756040961038 + 2.77555756156289e17*I, 0.124756040961038  1.38777878078145e17*I, 0.124756040961038 + 1.38777878078145e17*I, 0.124756040961038  2.77555756156289e17*I, 0.124756040961038 + 2.77555756156289e17*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.
comment:21 Changed 8 years ago by
 Commit changed from 3d9a07cd45fd002020a1f67cb2394af5a38d76a0 to 8274e4764945e773078535aabe9f262462416e57
Branch pushed to git repo; I updated commit sha1. New commits:
8274e47  Merge branch 'develop' into ticket/12947numerical_integral

comment:22 Changed 8 years ago by
 Status changed from needs_review to positive_review
As already said ...
comment:23 Changed 8 years ago by
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 1e16 Expected: 0.124756040961038 Got: 0.124756040961038  1.38777878078145e17*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 1e16 Expected: 0.124756040961038 Got: 0.124756040961038  2.77555756156289e17*I **********************************************************************
comment:24 Changed 8 years ago by
 Status changed from positive_review to needs_work
comment:25 Changed 8 years ago by
 Commit changed from 8274e4764945e773078535aabe9f262462416e57 to 72706eb6697fba361e138ddaca0765f8415eb8f3
Branch pushed to git repo; I updated commit sha1. New commits:
72706eb  Trac 12947: handle numerical noise correctly

comment:26 Changed 8 years ago by
 Status changed from needs_work to positive_review
This should solve it; probably no reason to ask for another review.
comment:27 Changed 8 years ago by
 Branch changed from u/pbruin/12947maxima_integral to 72706eb6697fba361e138ddaca0765f8415eb8f3
 Resolution set to fixed
 Status changed from positive_review to closed
To compare the "right" answer, do
Note how they don't really look quite similar, not really a branch cut thing? Also, in Maxima